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Abstract 

We review recent results on quantum reactive scattering from a phase space per- 
spective. The approach uses classical and quantum versions of Poincare-Birkhoff nor- 
mal form theory and the perspective of dynamical systems theory. Over the past ten 
years the classical normal form theory has provided a method for realizing the phase 
space structures that are responsible for determining reactions in high dimensional 
Hamiltonian systems. This has led to the understanding that a new (to reaction 
dynamics) type of phase space structure, a normally hyperbolic invariant manifold 
(or, NHIM) is the "anchor" on which the phase space structures governing reaction 
dynamics are built, e.g. it is the classical analogue of the chemists notion of the 
"activated complex" and it is essential for the construction of a surface that divides 
reactants from products which has the "no-recrossing" property for trajectories and 
minimal flux. The quantum normal form theory provides a method for quantizing 
these phase space structures through the use of the Weyl quantization procedure. 
We show that this approach provides a solution of the time-independent Schrodinger 
equation leading to a (local) S-matrix in a neighborhood of the saddle point gov- 
erning the reaction. These results can be obtained for any dimensional system for 
which an accurate normal form can be computed, and it does not require numerical 
solution of the Schrodinger equation or the generation of any classical trajectories. 
It follows easily that the quantization of the directional flux through the dividing 
surface with the properties noted above is a flux operator that can be expressed in 
a "closed form". Moreover, from the local S-matrix we easily obtain an expression 
for the cumulative reaction probability (CRP), which is the essential ingredient for 
the computation of microcanonical reaction rates and thermal reaction rates. Signif- 
icantly, the expression for the CRP can be evaluated without the need to compute 
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classical trajectories. This is a by product of the quantization of classical phase space 
structures that govern "exact" classical dynamics. The quantization of the NHIM 
is shown to lead to the activated complex, and the lifetimes of quantum states ini- 
tialized on the NHIM correspond to the Gamov-Siegert resonances. We apply these 
results to the collinear nitrogen exchange reaction and a three degree-of-freedom sys- 
tem corresponding to an Eckart barrier coupled to two Morse oscillators. We end 
by describing some further challenges that are topics of current research, but where 
some preliminary results are known: corner-cutting tunneling, state-to-state reac- 
tion rates, the flux-flux autocorrelation function formalism and the convergence of 
the quantum normal form. We emphasize that this dynamical systems, phase space 
approach to quantum reactive scattering through the quantum normal form provides 
a completely new approach to the computation of the relevant quantum scattering 
quantities (e.g. CRP, resonances) which shows promise in leading to computationally 
efficient methods for "high dimensional" systems. 
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1 Introduction 



Over the past fifty years the computation of quantities describing quantum reactive 
scattering processes has been a topic of great interest in the chemistry and physics 
community, see the Perspective article of Miller ^ for background and history. The 
interest in this topic continues to grow as a result of the need to include quantum 
effects in order to understand the behavior of certain biomolecules and molecular ma- 
terials (PJ). These contemporary applications are concerned with "large systems", 
and quantum mechanical calculations are notorious for their consumption of compu- 
tational resources as the number of degrees-of-freedom grows. As noted in |3j, the 
Nobel Prize winner Walter Kohn has identified the problem as the exponential wall of 
difficulty when one tries to perform numerical calculations in the traditional manner 
using the rules of quantum mechanics. In this review we discuss a new (phase space) 
approach to quantum reactive scattering that allows us to bypass the traditional 
computational difficulties for certain types of physically relevant quantities. 

Our purpose here is not to review the entire field of quantum reactive scattering. 
That would be a huge task that is beyond the scope of this review. Rather, we 
begin by describing the standard quantities that are computed to describe quantum 
reactive scattering processes and the issues associated with their computation. Our 
description follows the excellent review of Miller ^ . 

The quantity that contains a great deal of detailed information about a particular 
quantum reactive scattering process is the S-matrix. The S-matrix, {S'^^^n,, (S, J)} 
as a function of total energy E and total angular momentum J, can be obtained by 
solving the Schrodinger equation with scattering boundary conditions. From the S- 
matrix state-to-state differential scattering cross-sections can be computed as follows: 



where (rip) labels the reactant (product) rotational and vibrational states, 9 is the 
scattering angle between the relative velocity vectors of reactants and products, 
(mp) is the projection of total angular momentum onto the relative velocity vector 
of the reactants (products) and d:^^,{9) is the Wigner rotation matrix. 

For many chemical reactions a sufficient description is provided by the rate con- 
stant, either the canonical rate constant characterized by the temperature, k{T), or 
the microcanonical rate constant characterized by the total energy, k{E). These rate 
constants can be obtained using appropriate averages of the state-to-state differen- 
tial cross sections. This averaging process yields the cumulative reaction probability 



(1.1) 



J 



(CRP): 




(1.2) 
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which is used to compute the micro canonical and canonical rate constants, as follows 
k{E) = [2TThpr{E)]-^N{E), (1.3) 

/oo 
dEe-^^M{E). (1.4) 
-oo 

where /3 = {kT)~^ , is the density of reactant states per unit energy, and Qr is the 
reactant partition function per unit volume. 

Miller has pointed out ([4]) that if one is interested only in the rate, then solving 
the Schrodinger equation to obtain the S-matrix from which the state-to-state dif- 
ferential scattering cross section (jl.ip is computed, and then subsequently averaged 
to obtain the cumulative reaction probability (|1.2p . is an extremely inefficient pro- 
cess in the sense that a great deal of the detailed information obtained through the 
computation of the S-matrix is "thrown away" in the averaging process. Instead, 
Miller [3] discusses a direct way to compute M{E), or k(T), i.e. a method that avoids 
first computing the S-matrix, that is also correct in the sense that the expressions 
for J\f{E), or k{T) are exact, i.e. no approximations are involved. These expressions 
are given as follows (t4j): 

M{E) = 2TThTr(^5{E-H)FPr^ (1.5) 

k{T) = Qr{T)-'TT [e~^"FPr) (1.6) 

where H is the Hamilton operator, F is a flux operator, is the long time limit 
of the quantum mechanically time evolved Heaviside function, and Tr(-) denotes the 
trace operation. Schatz and Ratner [5] point out that there are three possibilities 
for the use of (jl.Sp and ()1.6p for computing rates. The ffi'st is to evaluate the trace 
of the operators in a basis, the second is to develop semiclassical theories for the 
rate constants ( [H [2] ) , and the third is to use the expressions to develop a version 
of quantum transition state theory ([7J). An alternative approach is to express p.Sp 
and ()1.6p in terms of a form of flux- flux autocorrelation function ([319]), and we 
discuss this approach in some detail in Section 16.31 A recent survey of methods for 
computing reaction rates is given in [lOj . 

Our approach us firmly rooted in phase space, both the classical and quantum 
theories. The classical phase space theory forms the "skeleton" on which our quan- 
tum theory of reactive scattering in phase space is built. The classical theory of 
reaction dynamics in phase space is reviewed in Section [2j Certainly there has been 
earlier work on classical reaction dynamics in phase space. In particular, we note the 
seminal work of Pechukas, Pollak, Child, and McLafferty from the 1970's and early 
1980's on two degree-of-freedom Hamiltonian systems f pTl [T2l [T3l [Tl [T5l [T6l [T71 118] ^ 
where the notion of a periodic orbit dividing surface (PODS) was introduced. The 
PODS bound a two dimensional dividing surface (in the three dimensional energy 
surface) having the "no-recrossing" property and minimal flux. This work did not 
generalize to more than two degrees-of-freedom. The generalization has required a 
new concept-the notion of a normally hyperbolic invariant manifold, or NHIM. This 
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was introduced in the study of phase space transport in Hamiltonian systems with 
three or more degrees-of- freedom in [19], and was subsequently used in the study of 
chemical reaction dynamics in |20] . For a d degree-of-freedom Hamiltonian system, 
the NHIM has the structure oi a 2d — 3 dimensional sphere in the 2d — 1 dimen- 
sional energy surface (so for d = 2, the NHIM has the structure of a one dimensional 
sphere, or periodic orbit). Conceptually, the NHIM is a fundamentally new dynami- 
cal object. It is a manifold of "saddle type stability" that exists in phase space. In 
Section [2] we discuss in more detail why this is the mathematical manifestation of the 
chemists notion of the activated complex. Later, it was shown in [21] that the NHIM 
played the role of the PODS for systems with three or more degrees-of-freedom in the 
sense that the NHIM is the boundary of a dividing surface having the "no-recrossing" 
property and minimal flux. 

This approach to phase space reaction dynamics for three or more degree-of- 
freedom Hamiltonian systems did not receive a great deal of applications throughout 
the 1990's. One reason for this was that there was not a computational method 
for realizing NHIMs in specific Hamiltonian systems (with three or more degrees-of- 
freedom) . This changed as a result of the work in [22\ [23] where it was shown that 
the classical Poincare-Birkhoff normal form theory could be used to compute a coor- 
dinate system where the NHIM could be realized, along with a variety of associated 
phase space structures. Software was developed that enabled the computation of the 
classical and quantum normal forms, and this allowed the realization of these phase 
space structures in specific exampleqj. This enabled the application of this approach 
to classical phase space reaction dynamics to HON isomerization ([23]), computation 
of the volume of points in an energy surface corresponding to reactive trajectories 
(the "reactive volume", |25l [26]), and a realization of Thiele's [27] notion of "gap 
times" and reaction rates ([28j). 

The generalization of this approach to the quantum setting has been carried 
out in [291 Sni [21] • The essential tools that make this possible are a quantum me- 
chanical generalization of the classical Poincare-Birkhoff normal form (described in 
Appendix lA.2p and the Weyl quantization procedure. This firmly establishes a phase 
space approach to quantum mechanics. Our development and understanding of clas- 
sical phase space structures associated with reaction dynamics is essential because 
one sees precisely their quantum mechanical manifestations. For example, in the 
classical setting the normal form is integrable in a (phase space) neighborhood of a 
particular saddle point that is relevant to the reaction of interest (in a way that we 
precisely describe in Section [2]). This property is manifested quantum mechanically 
by yielding a Hamilton operator in (quantum) normal form that factors into separate 
one-dimensional problems that can be solved exactly. This, in turn, gives rise to a 
block diagonal S-matrix that can easily be computed and analyzed. This is described 
in some detail in Section 14. li This is a significant result because the implication is 
if the relevant quantum normal form for a system can be computed, then regardless 
of the size of the system, the S-matrix can be calculated. In other words, our ap- 
proach leads to an expression for the S-matrix that does not require a solution of the 



^The software, along with documentation, is freely available for download at 
http : / /lacms . maths . bris . ac . uk/publications/ software/ index . html. 
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Schrodinger equation. In Section 14.11 we show that the classical flux through a di- 
viding surface in phase space having the no-recrossing and minimal flux properties is 
directly related to the quantum mechanical cumulative reaction probability. We ap- 
ply this approach to computing the cumulative reaction probability for the collinear 
nitrogen exchange reaction in Section 14.21 and coupled Eckart-Morse-Morse oscilla- 
tors in Section 14. 3[ In Section we show that quantum resonances are related to the 
time of decay of quantum states initialized on the activated complex. We conclude in 
Section E] by considering four developing areas of the theory: "corner cutting" tunnel- 
ing trajectories, state-to-state reaction rates, the flux-flux autocorrelation formalism, 
and convergence of the quantum normal form. 



2 Phase-Space Structures Underlying Reaction 
Dynamics 

In this section we describe the dynamics and geometry associated with phase space 
structures governing reaction dynamics. While the emphasis in this review is on 
quantum mechanics, the phase space structure that we describe forms the classical 
mechanical "skeleton" on which the quantum mechanical theory is supported. Here 
we merely summarize the basic results, more details can be found in the references 
given in the introduction. Our exposition here follows |32j . 

We consider a Hamiltonian system with d degrees of freedom, phase space coor- 
dinates {q,p) G E.'^ X K,'^ and Hamiltonian function H. We assume that {qo,Po) is an 
equilibrium point of Hamilton's equations which is of saddle-centre-. . .-centre stabil- 
ity type. El By adding a constant term to the Hamiltonian function (which does not 
change the dynamics) we can, without loss of generality, assume that H{qQ,pQ) = 0. 
Moreover, for simplicity of exposition, we can assume that the coordinates have been 
suitably translated so that the relevant equilibrium point {qo,Po) is at the origin. For 
much of the discussion below, we will consider iso-energetic geometrical structures 
belonging to a single positive energy surface := H~^{E) for some constant 

E > However, we note that in [33j essentially the same "picture" of phase space 
structures governing reaction dynamics has been shown to occur for constant tem- 
perature dynamics governed by the Hamiltonian isokinetic thermostat. 

We construct geometric structures in a neighborhood of the equilibrium point. 
We emphasize that by "neighborhood" we mean a neighborhood in phase space, 
not just on a fixed energy surface. We denote this neighborhood of the equilibrium 
point by C. We will defer until later a discussion of the constraints on the size of 
this neighborhood; suffice it to say for now that the region is chosen so that a new 
set of coordinates can be constructed (the normal form coordinates) in which an 
accurate nonlinear approximation of the Hamiltonian can be expressed (the normal 

^We will define this more precisely shortly. However, briefly, it nieans that the matrix associated with 
the linearization of Hamilton's equations about this equilibrium point has two real eigenvalues of equal 
magnitude, with one positive and one negative, and d — 1 purely imaginary complex conjugate pairs of 
eigenvalues. We will assume that the eigenvalues satisfy a nonresonance condition that we will describe 
more fully in the following. 
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form Hamiltonian) such that it provides an integrable approximation to the dynamics, 
as well as an algorithmic procedure for realizing phase space structures to within a 
given desired accuracy. 

Before describing the dynamics in normal form coordinates, as well as the realiza- 
tion of the phase space structures in normal form coordinates that govern "reaction" , 
i.e. trajectories in a phase space neighborhood of the saddle-centre-. . .-centre equilib- 
rium point (henceforth, referred to as a "saddle") we describe the relevant geometric 
structures that have been developed in the in the references given above. We em- 
phasize that although the normal form procedure provides a specific method for 
realizing these geometric structures in a given coordinate system, their existence is 
independent of any specific set of coordinates. 

Locally, the {2d— l)-dimensional energy surface ^{E) has the structure of S^*^^^ x 
K, in the 2(i-dimensional phase space. The energy surface is split locally into 

two components, "reactants" and "products" , by a {2d — 2)-dimensional "dividing 
surface" that is diffeomorphic to S^*^"^ and which we therefore denote by S'^^~'^{E). 
The dividing surface that we construct has the following properties:- 

• The only way that trajectories can evolve from reactants to products (and vice- 
versa), without leaving the local region C, is through S'^f~'^{E). In other words, 
initial conditions (ICs) on this dividing surface specify all reacting trajectories. 

• The dividing surface that we construct is free of local recrossings; any trajectory 
which crosses it must leave the neighbourhood £ before it might possibly cross 
again. 

• The dividing surface that we construct minimizes the flux, i.e. the directional 
flux through the dividing surface will increase upon a generic deformation of 
the dividing surface (see [21] for the details). 

The fundamental phase space building block that allows the construction of a 
dividing surface with these properties is a particular Normally Hyperbolic Invariant 
Manifold (NHIM) which, for a fixed positive energy E, win be denoted S^'^^^l^{E). 
The NHIM is diffeomorphic to S^'^"^ and forms the natural dynamical equator of 
the dividing surface: The dividing surface is split by this equator into {2d — 2)- 
dimensional hemispheres, each diffeomorphic to the open {2d — 2)-ball, B'^'^~^ . We 
will denote these hemispheres by B'^~^^{E) and B'^~^{E) and call them the "forward 
reactive" and "backward reactive" hemispheres, respectively. B'^~^^{E) is crossed 
by trajectories representing "forward" reactions (from reactants to products), while 
B'^~t^{E) is crossed by trajectories representing "backward" reactions (from products 
to reactants). 

The {2d — 3)-dimensional NHIM can be viewed as the energy surface of an (un- 
stable) invariant subsystem which as mentioned above, in chemistry terminology, 
corresponds to the "activated complex", which as an oscillating "super molecule" is 
located between reactants and products. 

The NHIM is of saddle stability type, having {2d — 2)-dimensional stable and 
unstable manifolds W'^{E) and W^{E) that are diffeomorphic to S^*^"^ x M. Being of 



8 



co-dimension □ one with respect to the energy surface, these invariant manifolds act 
as separatrices, partitioning the energy surface into "reacting" and "non-reacting" 
parts as will explain in detail in Sec. I2.2[ 

2.1 Dynamics and phase space structures near the sad- 
dle in normal form coordinates 

As mentioned in the previous section, reaction type dynamics are induced by equi- 
librium points of saddle X centre X . . . x centre stability type. These are equilibria for 
which the matrix associated with the linearisation of Hamilton's equations have eigen- 
values which consist of a pair of real eigenvalues of equal magnitude and opposite 
sign, (+A,— A), A G IR, and (d — 1) pairs of complex conjugate purely imaginary 
eigenvalues, —iuk), £ ioi k = 2, . . . , d. 

The phase space structures near equilibria of this type exist independently of a 
specific coordinate system. However, in order to carry out specific calculations we 
will need to be able to express these phase space structures in coordinates. This is 
where Poincare-Birkhoff normal form theory is used. This is a well-known theory and 
has been the subject of many review papers and books, see, e.g., [Ml ESI 133 EZl ESI 
[391 Bo] . For our purposes it provides an algorithm whereby the phase space structures 
described in the previous section can be realised for a particular system by means of 
the normal form transformation which involves making a nonlinear symplectic change 
of variables. 



into normal form coordinates, {q, p) which, in a local neighbourhood C of the equi- 
librium point, "decouples" the dynamics into a "reaction coordinate" and "bath 
modes." The coordinate transformation T is obtained from imposing conditions on 
the form of H expressed the new coordinates, {q,p), 



These conditions are chosen such that -ffcNF and the resulting equations of motions 
assume a simple form in which the reaction coordinate and bath modes "decouple" . 
This decoupling is one way of understanding how we are able to construct the phase 
space structures, in the normal form coordinates, that govern the dynamics of reac- 
tion. 

In fact, we will assume that a (generic) non-resonance condition holds between 
the eigenvalues, namely that 



■^Briefly, the co-dimension of a submanifold is the dimension of the space in which the submanifold 
exists, minus the dimension of the submanifold. The significance of a submanifold being "co-dimension 
one" is that it is one less dimension than the space in which it exists. Therefore it can "divide" the space 
and act as a separatrix, or barrier, to transport. 



{q,p) = T{q,p) 



(2.1) 




(2.2) 



k2UJ2 H h knUJd / 



(2.3) 
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for all integer vectors {k2, ■ ■ ■ ,kd) S Z'^^^. When such a condition holds, the 
classical normal form procedure, see Appendix lA.H yields an explicit expression for 
the normalised Hamiltonian Hqtssf as a function of d integrals of motion: 

d [_N/2\ 

= Eo + XI + ^ujkJk+ ^ Yl -^2' ■ ■ ■ Jd' ■ (2-4) 

k=2 n=2 \a\=n 

Note that the normal form is expressed as a polynomial in the action integrals / and 
Jfc, k = 2, . . . ,d, corresponding to a "reaction coordinate" and {d — 1) "bath modes" 
respectively. In Eq. (|2.4p £'0 denotes the energy of the system at the equilibrium 
point, and the expansion coefficient obtained by the classical normal form 

algorithm presented in Appendix [Xj here, a = (ai, . . . ,0^) with G Mq for k = 
1, . . . ,d, and |a| = a^. In Appendix|A] we show that the normal form algorithm 
generates a polynomial at each step, where the normalization algorithm at a certain 
order does not modify the normalized terms at lower order. However, in practice it 
is necessary to stop the algorithm at a given order. In this sense the normal form is 
an approximation to the original Hamiltonian in a neighborhood of the saddle, which 
we have denoted by C. Certainly, it is essential that the approximation is accurate 
enough to yield useful information. Numerous examples have show this to be the 
case, and we discuss this issue more fully in Section 12.51 

The integral, /, corresponds to a "reaction coordinate" (saddle- type DoF): 

I = qiPl. (2.5) 

We note that there is an equivalent form of the reaction coordinate:- making the 
linear symplectic change of variables qi = (Pi — Qi)/^/2 and pi = {Pi + Qi)/\/2, 
transforms the above into the following form, which may be more familiar to many 
readers, 

I = qm = l{P?-Ql). (2.6) 

Geometrically speaking, one can move freely between these two representations by 
considering the plane {qi,pi) and rotating it by angle 7r/4, to give {Q,P). 

The integrals J^, for k = 2, . . . ,d, correspond to "bath modes" (centre-type DoF 

s 

"'We note that the inclusion of ±A in a non-resonance condition would be vacuous; one cannot have a 
resonance of this kind between a real eigenvalue, ±A, and purely imaginary eigenvalues, izwfe, A: = 2, . . . , d. 

'"^Throughout our work we use, somewhat interchangeably, terminology from both chemical reaction 
dynamics and dynamical systems theory. This is most noticable in our reference to the integrals of 
motion. / is the integral related to reaction, and in the context of dynamical systems theory it is related 
to hyperbolic behaviour. The term "reactive mode" might also be used to describe the dynamics associated 
with this integral. The integrals J2, ■ ■ ■ , Jd describe the dynamics associated with "bath modes". In the 
context of dynamical systems theory, the dynamics associated with these integrals is referred to as "center 
type dynamics" or "center modes" . A key point here is that integrals of the motion provide us with the 
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Jk = l{pl + ql), k = 2,...,d. (2.7) 
In the new coordinates, Hamilton's equations have a particularly simple form: 



dHcNF \ IT 1 T \ 



dp 



Pi = 7\ = -HI, J2,---, Jd)Pl, 



dq 



1 



qk — — ^ — "fciJ ) <^2, • • • , -'djPfc 



(2. 



k 



Pk — ^ — -"fci-'; ^2, ■ ■ ■ , Jdjqk, 



dqk 

for k = 2, . . . ,d, where we denote 



A(/,J„...,J.) := 9KcMl,J2,...,Jd) (2.9) 
n,iI,J,,...,J,) := 9KcMY2,...,J,) ^ k = 2,...,d. (2.10) 

The integrals provide a natural definition of the term "mode" that is appropriate 
in the context of reaction, and they are a consequence of the (local) integrability 
in a neighborhood of the equilibrium point of saddle-centre-. . .-centre stability type. 
Moreover, the expression of the normal form Hamiltonian in terms of the integrals 
provides us a way to partition the "energy" between the different modes. We will 
provide examples of how this can be done in the following. 

The normal form transformation T in ()2.ip can be computed in an algorithmic 
fashion. One can give explicit expression for the phase space structures discussed in 
the previous section in terms of the normal form coordinates, {q,p)- This way the 
phase space structures can be constructed in terms on the normal form coordinates, 
{q,p), and for physical interpretation, transformed back to the original "physical" 
coordinates, {q,p), by the inverse of the transformation T. H 

In summary, the "output" of the normal form algorithm is the following: 

• A symplectic transformation T{q,p) = {q,p), and its inverse T~^{q,p) = {q,p), 
that relate the normal form coordinates {q,p) to the original "physical" coor- 
dinates {q,p)- 



natural way of defining and describing the physical notion of a "mode" . The nature of the mode is defined 
in the context of the specific application (i.e. chemical reactions) or, in the context of dynamical systems 
theory, through its stability properties (i.e. hyperbolic or centre). 

®The original coordinates {q,p) typically have an interpretation as configuration space coordinates and 
momentum coordinates. The normal form coordinates {q,p), in general, do not have such a physical 
interpretation since both q and p are nonlinear functions of both q and p. 
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• An expression for the normalized Hamiltonian: in the form, f^cNFj in terms 
of the normal form coordinates {q,p), and in the form JCcNFj in terms of the 
integrals (/, J2, . . . , Jd)- 

• Explicit expressions for the integrals of motion / and J^, k = 2, . . . ,d, in terms 
of the original "physical" coordinates by means of the symplectic transformation 
T{q,p) = {q,p). 

2.2 Explicit definition and construction of the phase space 
structures in the normal form coordinates 

As indicated in the previous section it is straightforward to construct the local phase 
space objects governing "reaction" in the normal form coordinates, {q,p). In this 
section, we will define the various structures in the normal form coordinates and 
discuss briefly the consequences for the original dynamical system. 

The structure of an energy surface near a saddle point: For E < 0, the 

energy surface consists of two disjoint components. The two components correspond 
to "reactants" and "products." The top panel of Fig. [T] shows how the two compo- 
nents project to the various planes of the normal form coordinates. The projection 
to the plane of the saddle coordinates {qi,pi) is bounded away from the origin by 
the two branches of the hyperbola, qipi = / < 0, where / is given implicitly by 
the energy equation with the centre actions J^, k = 2,...,d, set equal to zero: 
-f^CNF(-^;0, . . . ,0) = E < 0. The projections to the planes of the centre coordinates, 
{qk,Pk), k = 2, . . . ,d, are unbounded. 

At -E = 0, the formerly disconnected components merge (the energy surface bi- 
furcates), and for E > 0, the energy surface has locally the structure of a spherical 
cylinder, S^*^"^ x IR,. Its projection to the plane of the saddle coordinates now includes 
the origin. In the first and third quadrants it is bounded by the two branches of the 
hyperbola, qipi = / > 0, where / is again given implicitly by the energy equation 
with all centre actions equal to zero, but now with an energy greater than zero: 
-f^CNF(-^, 0, . . . ,0) = E > 0. The projections to the planes of the centre coordinates 
are again unbounded. This is illustrated in the bottom panel of Fig. [TJ 

The dividing surface, and reacting and nonreacting trajectories: On an 

energy surface with £^ > 0, we define the dividing surface by qi = pi- This gives a 
(2(i— 2)-sphere which we denote by S'^'^~'^{E). Its projection to the saddle coordinates 
simply gives a line segment through the origin which joins the boundaries of the 
projection of the energy surface, as shown in Fig. [21 The projections of the dividing 
surface to the planes of the centre coordinates are bounded by circles {p1+q'^)/2 = J^, 
k = 2, . . . ,d, where is determined by the energy equation with the other centre 
actions, Ji, I ^ k, and the saddle integral, /, set equal to zero. The dividing surface 
divides the energy surface into two halves, pi — qi > and pi — qi < 0, corresponding 
to reactants and products. 

As mentioned above, trajectories project to hyperbolae in the plane of the saddle 
coordinates, and to circles in the planes of the centre coordinates. The sign of / 
determines whether a trajectory is nonreacting or reacting, see Fig. [2j Trajectories 
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Figure 1: Projection of energy surfaces (turquoise regions) to the planes of the normal form 
coordinates. The energy surface have energy E < (top panel), E = (middle panel), 
E > (bottom panel). 



which have / < are nonreactive and for one branch of the hyperbola qipi = I they 
stay on the reactants side and for the other branch they stay on the products side; 
trajectories with / > are reactive, and for one branch of the hyperbola qipi = I 
they react in the forward direction, i.e., from reactants to products, and for the other 
branch they react in the backward direction, i.e., from products to reactants. The 
projections of reactive trajectories to the planes of the centre coordinates are always 
contained in the projections of the dividing surface. In this, and other ways, the 
geometry of the reaction is highly constrained. There is no analogous restriction on 
the projections of nonreactive trajectories to the centre coordinates. 

The normally hyperbolic invariant manifold (NHIM) and its relation to 
the "activated complex": On an energy surface with E > 0, the NHIM is given 
by Qi = Pi = 0. The NHIM has the structure of a (2c? — 3)-sphere, which we denote 
by 5'n'^iim(^)- The NHIM is the equator of the dividing surface; it divides it into 
two "hemispheres": the forward dividing surface, which has qi = pi > 0, and the 
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P2, 



P3 



Figure 2: Projection of the dividing surface and reacting and nonreacting trajectories to 
the planes of the normal form coordinates. In the plane of the saddle coordinates, the 
projection of the dividing surface is the dark red diagonal line segment, which has qi = pi- 
In the planes of the centre coordinates, the projections of the dividing surface are the 
dark red discs. Forward and backward reactive trajectories (yellow and blue) project to 
the first and third quadrant in the plane of the saddle coordinates, respectively, and pass 
through the dividing surface. The red and green curves mark nonreactive trajectories on 
the reactant side (pi — qi > 0), and on the product side {pi — gi < 0), of the dividing 
surface, respectively. The turquoise regions indicate the projections of the energy surface. 



backward dividing surface, which has qi = pi < 0. The forward and backward 
dividing surfaces have the structure of {2d — 2)-dimensional balls, which we denote 
by B'^~^^{E) and B'^f~^{E), respectively. All forward reactive trajectories cross 

-^ds7^(-^)' backward reactive trajectories cross B^~t^{E). Since gi = pi = in 
the equations of motion (|2.8p implies that qi = "Pi = 0, the NHIM is an invariant 
manifold, i.e., trajectories started in the NHIM stay in the NHIM for all time. The 
system resulting from q\ = pi = Q \s an invariant subsystem with one degree of 
freedom less than the full system. In fact, qi = pi = Q defines the centre manifold 
associated with the saddle-centre-- • • -centre equilibrium point, and the NHIM at an 
energy E greater than the energy of the quilibrium point is given by the intersection 
of the centre manifold with the energy surface of this energy E |231 [2T] . 

This invariant subsystem with one degree of freedom less than the full system is 
the "activated complex" (in phase space), located between reactants and products. 
The NHIM can be considered to be the energy surface of the activated complex. In 
particular, all trajectories in the NHIM have / = 0. 

The equations of motion (j2.8p also show that pi — qi < on the forward divid- 
ing surface B'^~^^{E), and pi — gi > on the backward dividing surface B'^~^{E). 
Hence, except for the NHIM, which is is an invariant manifold, the dividing surface is 
everywhere transverse to the Hamiltonian flow. This means that a trajectory, after 
having crossed the forward or backward dividing surface, B'^~^{E) or B'^~^{E), 
respectively, must leave the neighbourhood of the dividing surface before it can pos- 
sibly cross it again. Indeed, such a trajectory must leave the local region in which 
the normal form is valid before it can possibly cross the dividing surface again. 

The NHIM has a special structure: due to the conservation of the centre actions. 
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it is filled, or foliated, by invariant {d — l)-dimensional tori, T'^^^. More precisely, 
for d = 3 degrees of freedom, each value of J2 implicitly defines a value of J3 by the 
energy equation i^cNF(0, J2, -^s) = E. For three degrees of freedom, the NHIM is 
thus foliated by a one-parameter family of invariant 2-tori. The end points of the 
parameterization interval correspond to J2 = (implying ^2 = P2 = 0) and J3 = 
(implying qs = P3 = 0), respectively. At the end points, the 2-tori thus degenerate 
to periodic orbits, the so-called Lyapunov periodic orbits. 




Figure 3: The projection of the NHIM and the local parts of its stable and unstable 
manifolds, W^{E) and W^lE), to the planes of the normal form coordinates. In the plane 
of the saddle coordinates, the projection of the NHIM is the origin marked by the blue bold 
point, and the projection of W^{E) and W'^{E) are the pi-axis and gi-axis, respectively. 
W^{E) consists of the forward and backward branches Wf{E) and W^{E), which have 
Pi > and pi < 0, respectively; W^{E) consists of Wl^{E) and W^{E), which have gi > 
and qi < 0, respectively. In the plane of the centre coordinates, the projections of the 
NHIM, W^{E), and W^{E) (the blue circular discs) coincide with the projection of the 
dividing surface in Fig. [2l The turquoise regions mark the projections of the energy surface. 



The stable and unstable manifolds of the NHIM forming the phase space 
conduits for reactions: Since the NHIM is of saddle stability type, it has stable 
and unstable manifolds, W^{E) and W^{E). The stable and unstable manifolds have 
the structure of spherical cylinders, S^*^"^ x R. Each of them consists of two branches: 
the "forward branches", which we denote by Wj{E) and W^{E), and the "backward 
branches", which we denote by W^{E) and W^{E). In terms of the normal form 
coordinates, Wl{E) is given by gi = with pi > 0, W^{E) is given by = with 
qi > 0, W^{E) is given by (?! = with pi < 0, and W^^^E) is given by pi = with 
qi < 0, see Fig. [3l Trajectories on these manifolds have / = 0. 

Since the stable and unstable manifolds of the NHIM are of one less dimension 
than the energy surface, they enclose volumes of the energy surface. We call the 
union of the forward branches, Wj{E) and W^{E), the forward reactive spherical 
cylinder and denote it by Wi{E) . Similarly, we define the backward reactive spherical 
cylinder, Wb{E), as the union of the backward branches, W^{E) and W^{E). 

The reactive volumes enclosed by Wf{E) and Wb{E) are shown in Fig. [5] as 
their projections to the normal form coordinate planes. In the plane of the saddle 
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Figure 4: Projections of the reactive volumes enclosed by the forward and backward reac- 
tive spherical cylinders, Wf{E) and Wh{E), and the forward and backward reactions paths, 
to the planes of the normal form coordinates. The volumes enclosed by W{{E) and Wh{E) 
project to the dark pink and green regions in the first and third quadrant in the plane of 
the saddle coordinates, respectively. These volumes project to the dark green/dark pink 
brindled disks in the planes of the centre coordinates, where their projections coincide with 
the projection of the NHIM and the dividing surface in Figs. [2] and [31 The forward and 
backward reaction paths project to the two branches of a hyperbola marked blue in the 
first and third quadrant in the plane of the saddle coordinates, respectively, and to the 
origins (bold blue points) in the planes of the centre coordinates. The turquoise regions 
mark the projections of the energy surface. 

coordinates, the reactive volume enclosed by Wf{E) projects to the first quadrant. 
This projection is bounded by the corresponding hyperbola qipi = I, with I obtained 
from Kcnf(-^, 0, . . . , 0) = E. Likewise, Wb{E) projects to the third quadrant in 
the (gi,pi)-plane. Wf{E) encloses a// forward reactive trajectories; Wb{E) encloses 
all backward reactive trajectories. All nonreactive trajectories are contained in the 
complement. 

2.3 The foliation of the reaction region by Lagrangian 
submanifolds 

The existence of the d integrals of motion, (/, J2, . . . , Jrf), induce phase space struc- 
tures which lead to further constraints on the trajectories in addition to the ones 
described above. In order to describe these structures and the resulting constraints 
it is useful to introduce the so called momentum map, A4 ^41i. i42j which maps a point 
{qi, . . . ,qn,pi, ■ ■ ■ ,Pd) in the phase space R*^ x M!^ to the d integrals evaluated at this 
point: 

7W : X R'^ ^ {qi, . . . ,qn,pi, . . . ,Pd) ^ iI,J2, ■ ■ ■ ,Jd) ■ (2.11) 

The preimage of a value for the constants of motion (/, J2, . . . , Jd) under 7W is called 
a fibre. A fibre thus corresponds to the common level set of the integrals in phase 
space. 
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A point (gi, . . . , qd,Pi, ■ ■ ■ ,Pd) is called a regular point of the momentum map if 
the linearisation of the momentum map, DA^, has full rank d at this point, i.e., if 
the gradients of the n integrals /, J^, k = 2, . . . , d, with respect to the phase space 
coordinates {q,p) are linearly independent at this point. If the rank of DA^ is less 
than d then the point is called irregular. A fibre is called regular if it consists of 
regular points only. Otherwise, it is called an irregular fibre. In fact almost all fibres 
are regular. They are d-dimensional manifolds given by the Cartesian product of an 
hyperbola qipi = I in the saddle plane {qi,Pi) and d—1 circles in the centre planes 
{qk,Pk), k = 2, . . . ,d. Since the hyperbola qipi = I consists of two branches each of 
which have the topology of a line M, the regular fibres consist of two disjoint toroidal 
cylinders, T'^~^ x R, which are the Cartesian products of a (d — l)-dimensional torus 
and a line. We denote these toroidal cylinders by 



Kj2,-,J, = ^ ^'"^ • ^1^1 = ^' ^(^2 + Q2) = J2,--- ,lipl + ql) =Jd,qi>0} 

(2.12) 

and 



^7, J2,...,J, e Et'"' : Piqi=I, ^{Pl + q^) =J2,... ,^{pl + qd) = Jd,qi < 0} . 

(2.13) 

A^j^ Jd J2 Jd ^"^^ Lagrangian manifolds [43j. The Lagrangian manifolds 

consists of all trajectories which have the same constants of motion. In particular 
the Lagrangian manifolds are invariant, i.e. a trajectory with initial condition on a 
Lagrangian manifold will stay in the Lagrangian manifold for all time. For / < 0, the 
Lagrangian manifolds and A^j^ consist of nonreactive trajectories in 

the reactants resp. products components of the energy surface. For / > 0, A^ 
consists of forward reactive trajectories, and A^ consists of backward reactive 

trajectories. 



2.4 The Directional Flux Through the Dividing Surface 

A key ingredient of transition state theory and the classical reaction rate is the direc- 
tional flux through the dividing surface defined in Sec. 12. 2[ Given the Hamiltonian 
function in normal form expressed as a function of the integrals ()2.4p . and a fixed 
energy E above the energy of the saddle-centre-- •• -centre, Eq, it is shown in [21j 
that the directional flux through the dividing surface is given by 



f{E) = {2TTf-'V{E), (2.14) 

where V{E) is the volume in the space of the actions (J2, . . . , Jd) enclosed by the 
contour i^'cNF(0, J2, • • • , -/d) = E. This is a significant result because it enables 
the computation of the directional flux without computing trajectories. Moreover, it 
directly connects the directional flux to the NHIM, i.e. the activated complex. In 
Fig. [5] we illustrate the volume V{E) for the case of a 3 DoF systems: here V(-E') is 
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given by the area in the (J2, ./s) plane enclosed by the light blue line corresponding 
to the NHIM. 



J3A 




Figure 5: Contour i^'cNF(0, J2, ■ ■ ■ , Jd) = E (blue line) in the space of the centre integrals 
( J2, . . . , Jd) for (i = 3 degrees of freedom. Up to the prefactor {2nY~^, the area V{E) of 
the enclosed region (marked green) gives the directional flux through the dividing surface, 
see Equation (12.141) . The green region is the projection of the / > piece of the energy 
surface on the ( J2, J3)-plane. 



We also note here that the dimensionless quantity 



Mveyl(S)= (^^j^, (2.15) 



where 2TTh is Planck's constant, is Weyl's approximation of the integrated density of 
states, or equivalently the mean number of quantum states of the activated complex 
with energies less than or equal to E (see, e.g. [44J)- It is shown in [30J that A/'weyi(-E') 
can be interpreted as the mean number of open quantum "transition channels" at 
energy E. 

In the case where we only take into account the quadratic part of the normal form, 
or equivalently, if we linearise Hamilton's equations, we have Kc^p(I, J2, . . . , Jd) = 
XI + Ylk=2'^kJk and the energy surface -ftrcNF(0, J2, . . . , Jd) = E encloses a simplex 
in ( J2, . . . , Jd) whose volume leads to the well-known result 



This shows, e.g, that the flux scales with E'^~^ for energies close to the saddle energy. 
The key advantage of the normal form coordinates resulting from the normal form 
algorithm is that it allows one to include the non-linear corrections to ()2.16p to any 
desired order. 
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2.5 Effect of Truncation of the Normal Form Algorithm 

The normalization procedure proceeds via formal power series manipulations whose 
input is a Taylor expansion of the original Hamiltonian, H, necessarily up to some 
finite order, A^, in homogeneous polynomials. For a particular application, this pro- 
cedure naturally necessitates a suitable choice of the order, A^, for the normalization, 
after which one must make a restriction to some local region, £, about the equi- 
librium point in which the resulting computations achieve some desired accuracy. 
Hence, the accuracy of the normal form as a power series expansion truncated at 
order iV in a neighborhood C is determined by comparing the dynamics associated 
with the normal form with the dynamics of the original system. There are several 
independent tests that can be carried out to verify accuracy of the normal form. 
Straightforward tests that we use are the following: 

• Examine how the integrals associated with the normal form change on trajec- 
tories of the full Hamiltonian (the integrals will be constant on trajectories of 
the normal form). 

• Check invariance of the different invariant manifolds (i.e. the NHIM and its 
stable and unstable manifolds) with respect to trajectories of the full Hamilto- 
nian. 

Both of these tests will require us to use the transformations between the original 
coordinates and the normal form coordinates. Specific examples where A'^, C and 
accuracy of the normal forms are considered can be found in [JH HSJ [251 HZl ES]. 

2.6 "Globalizing" the Geometrical Structures in the Re- 
action Region 

As we have shown, the normal form transformation to normal form coordinates pro- 
vides a method for providing a complete understanding of the geometry of reaction 
dynamics in a neighbourhood C (in phase space) of the saddle-centre-- • • -centre equi- 
librium point of Hamilton's equations. By this, we mean that in the normal form 
coordinates we can give an explicit equation for the surfaces and, as a result of 
the "simple" structure of Hamilton's equations in the normal form coordinates, we 
can describe precisely the influence of these geometrical structures on trajectories of 
Hamilton's equations. In Tab. [I] we summarize the results obtained this far by pro- 
viding a list of the different surfaces that control the evolution of trajectories from 
reactants to products in the neighbourhood C in Fig. [H 

However, it must be kept in mind that all of these surfaces, and associated dy- 
namical phenomena, are only "locally valid" in the neighbourhood C . The next step 
is to understand their influence on the dynamics outside of C , i.e., their influence 
on the dynamics of reaction throughout phase space in the original coordinates (as 
opposed to the normal form coordinates) . In order to do this we will need the normal 
form transformation discussed in Appendix lAl to order A^ (where N is determined 
according to the desired accuracy following the discussion in Section 12. 5p . 
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Geometrical Structure 


Equation in Normal Form Coordinates 


dividing surface, Sl^'^{E) 


Qi = Pi 


forward reactive hemisphere, B^~^^{E) 


qi = Pi> 


backward reactive hemisphere, B'i'^~Jf(E) 

^ ' as, b v / 


Ql = Pl < 


NHIM, S-^'^^UE) 


qi=Pi = 


stable manifold of the NHIM, W%E) 


(7, = 0, »i 7^ 


unstable manifold of the NHIM, WiE) 


Pl = 0, gi ^ 


forward branch of ^^{E), Wf{E) 


qr, = 0, »i > 


backward branch of W^E), W^{E) 


gi = 0, Pl < 


forward branch of WiE), Wf{E) 


= 0, gi > 


backward branch of WiE), W^{E) 


Pl = 0, gi < 


forward reactive spherical cylinder 
Wf{E) = W'f{E) U Wf{E) 


Pigi = 0, pi,gi > 0, gi 7^ pi 


backward reactive spherical cylinder 
WbiE) = W,'iE)UWf^iE) 


Pigi = 0, pi,gi < 0, gi 7^ pi 



Table I: Table of phase space surfaces influencing reaction dynamics and their representa- 
tions in normal form coordinates on an energy surface of energy greater than the energy 
of the saddle equilibrium point. 

In Appendix [A] we discuss the necessary transformations of the original physical 
coordinates required to transform the Hamiltonian into normal form. In particular, 
we translate the saddle-centre-- • • -centre equilibrium point to the origin, we "sim- 
plify" the linear part of Hamilton's equations (what "simplify" precisely means is 
described in the Appendix), then we iteratively construct a sequence of nonlinear 
coordinate transformations that successively "simplify" the order 3, 4, . . ., N terms 
of the Hamiltonian according to the algorithm described in Appendix [Xj We can 
invert each of these transformations to return from the normal form coordinates to 
the physical coordinates. 

Computation of W^{E) and W^{E): Our approach to computing the stable 
and unstable manifolds of a NHIM is, in principle, the same as for computing the 
stable and unstable manifolds of a hyperbolic trajectory (however, the practical im- 
plementation of the algorithm in higher dimensions is a different matter and one that 
deserves much more investigation). 

We describe the computation of Wj{E) as follows. 

• In the normal form coordinates, choose a distribution of initial conditions on 
the NHIM and displace these initial conditions "slightly" in the direction of the 
forward branch of W^{E) (pi = 0, qi = e > 0, e "small"). 

• Map these initial conditions back into the physical coordinates using the inverse 
of the normal form transformation. 
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• Integrate the initial conditions forward in time using Hamilton's equations in 
the physical coordinates, for the desired length of time (typically determined 
by accuracy considerations) that will give the manifold of the desired "size". 
Since the initial conditions are in the unstable manifold they will leave the 
neighbourhood £ in which the normal form transformation is valid (which is 
why we integrate them in the original coordinates with respect to the original 
equations of motion). 

The backward branch of W^{E) can be computed in an analogous manner by dis- 
placing the initial conditions on the NHIM in the direction of the backward branch 
of WiE) (pi = 0, gi = e < 0, £ "small"). 

Computation of W^{E) and Wj{E): The forward and backward branches of 
W^{E) can be computed in an analogous fashion, except the initial conditions are 
integrated backward in time. 

Practical considerations: By their very definition, invariant manifolds consist 
of trajectories, and the common way of computing them, and visualizing them, that 
works well in low dimensions is to integrate a distribution of initial conditions lo- 
cated on the invariant manifold (hence, this illustrates the value of the normal form 
coordinates and transformation for locating appropriate initial conditions). In high 
dimensions there are numerical and algorithmic issues that have yet to be fully ad- 
dressed. How does one choose a mesh on a 2(i — 3 dimensional sphere? As this mesh 
evolves in time, how does one "refine" the mesh in such a way that the evolved mesh 
maintains the structure of the invariant manifold? 

Examples where this "globalization" of phase space structures in the reaction 
region has been carried out can be found in [25] [26l [2^ I46j . However, there is 
tremendous scope for future work in this direction, both from the point of view of 
mathematical and computational techniques, and applications to chemical dynamics. 

3 Quantum Normal Form Representation of 
the Activated Complex 

In this section we present an extension of the classical normal form theory to quantum 
mechanics. Leaving all technical details for Appendix |A] we only consider here key 
aspects and results of the theory. 

The quantum normal form (QNF) theorjl^ provides an explicit algorithmic proce- 
dure that allows one to approximate the Hamiltonian operator H, corresponding to a 
reactive system with a generally nonintegrable classical Hamiltonian, by a "simpler" 
Hamiltonian operator Hq^p, having an integrable classical counterpart. The operator 
-f^QNF takes the form of a power series expansion in reactive and nonreactive (bath- 
mode) action operators and the effective Planck's constant. The (generally infinite) 

^See Ref. [3D] for a comprehensive review of the mathematical theory, or consult Appendix |^ for a 
concise summary. 



power series can then be terminated at any desired order dictated by the desired 
accuracy of the approximation. This leads to an A^th order QNF approximation 
^QNF original Hamiltonian H: 

-^QNF = ^qnf(^' J2,Jz,---,Jd) 

d lN/2\ 

= Eo + \i + Y,^kJk+ ^n,aj''' JT ■ ■ ■ Jd'f^ee' . (3.1) 

fc=2 n=2 \a\+j=n 

This approximation holds locally, in the vicinity of a single equilibrium point of the 
Hamiltonian function of the corresponding d-dimensional classical system (see Ap- 
pendix |A] for details). It is this equilibrium point whose phase-space neighborhood 
is traversed by reactive trajectories on their way from the valley of reactants to the 
valley of products. The equilibrium point is considered to be of saddle-center-... - 
center stability type, meaning that the 2d x 2d matrix associated with the linearized 
Hamilton's equations of motion has two real eigenvalues, itA, and d—1 purely imag- 
inary complex conjugate pairs of eigenvalues, Huj). with k = 2, . . . ,d. The linear 
frequencies ujk are further assumed to be rationally independent, so that the condi- 
tion m2W2 + . . . + rUfiOJd = implies ra2 = . . . = rrifi = for all integers m2, . . . , m^. 
Other quantities entering Eq. ()3.ip include (i) the energy £"0 of the system at the 
equilibrium point, (ii) parameters Kn,ai,...,«d,j obtained from the QNF expansion pro- 
cedure, (iii) an effective (scaled) Planck's constant ^eff) which in atomic units is given 
by the inverse of the square root of the reduced mass of the system and plays the 
role of a "small parameter" , (iv) an action integral operator 

i = \{qiPi+Piqi) (3.2) 
corresponding to the reactive mode, and (v) action integral operators 

Jk = \{ql+Pl) . k = 2,...,d (3.3) 

of the bath modes. Here, the "position" and "momentum" operators, qk and pk 
respectively, satisfy the commutation relations 

[qk, qi] = [Pk,Pi] = , [qk,Pi] = iKs hi , (3.4) 

where k,l = 1, . . . ,d, and 6^1 is the Kronecker's delta. The operators qi, pi, . . ., ijd, 
Pd are given by the Weyl quantization of classical phase-space coordinates qi, pi, . . . , 
qa, Pd that in turn are obtained from the original "physical" phase-space coordinates 
by way of a nonlinear canonical transformation j30]. Thus, qi and pi correspond to 
the reaction coordinate, while qk and pk to the A;th bath mode with k = 2, . . . ,d. 

It sometimes proves convenient to work in a phase-space coordinate basis (Qi, 
Pi, q2, P2-, ■ ■ qd, Pd) (which we further refer as to "QP-basis") that is "rotated" at 
an angle 7r/4 with respect to the basis {qi, pi, q2, P2, ■ ■ qd, Pd) (further referred as 
to "qp-basis"), i.e.. 
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In the QP-basis, the action integral operator corresponding to the reaction coordinate 
reads 



1 



Pi - Ql 



(3.6) 



The main advantage of having the original Hamiltonian H approximated by a 
polynomial in the operators / and J^, A: = 2, . . . , d, is that the eigenstates of the QNF 
operator can be chosen to be simultaneously the eigenstates of the operators 

/ and jfc, whose spectral properties are well known. Indeed, 



with 
where 



^qnfI^' f^2, ■ • • , rid) = E\I, n2, . . . , n^) 

|/, 712, ... , nd) = \lpl) O IV'na) O • • • tX" IV'nd) , 



and 



,(7V) 



/ G E, 
nfc G INo , 

,Ks{nd + 1/2) 



(3.7) 
(3.8) 

(3.9a) 
(3.9b) 

(3.10) 



^ = ^QNH^'^cfT(n2 + l/2), 

And, after having explicitly described the eigenstates {ipj) and \ipnk) of the reaction 
and bath degrees of freedom respectively, one obtains a complete (approximate) 
eigensystem for the original reactive scattering problem. 



/ > 




/ < 




Figure 6: The blue lines show the Lagrangian manifolds Af'\ aJ'^^ Af''\ and A^"'"^ 
associated with the J-eigenstates {ip^j'"^^), lip^j'"^^), {ipf'''^^), and IV'/"'^^) respectively. The 
arrows indicate the classical Hamiltonian vector fields generated by / = piQi- The red 
thick line corresponds to the dividing surface s{qi,pi) = qi — pi- 

Two independent solution of Eq. ()3.9ap . in the gi-representation, are given by 



{o,r)\ 



{o,p)\ 



■ 0(^i)gi 



(3.11a) 
(3.11b) 
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where G is the Heaviside step function, the superscript "o" stands for "outgoing to" , 
and "r" and "p" stand for "reactants" and "products" respectively. The motivation 
for this notation is clear from viewing the solutions given by Eq. (j3.9p as Lagrangian 
states, i.e., by writing 

{q^^''"^) = A^r^^'Hgi) exp > (3-12) 

where the purely real amplitude and phase functions are given by 

Af''/'\qi) = Q{Tgi)\qi\-'^' , = nn\q,\ , (3.13) 

respectively. This way one can associate the one-dimensional Lagrangian manifolds 

A^''^ = [iqi,Pi) = (qi, -^/rHu)) = (ft, ^) : ft < o} , (3.14a) 

A?'^^ = [iqi,Pi) = (l^^^/rhQi)) = (ft>^) : ft > o} (3.14b) 

with the states I'lpf''^^) and 1-0/°'^^). From the presentation of A^^'^^ and Aj°'^^ in 
Fig. [6] one sees that for qi — )• — oo the wave function {qi\Tpf'''^^) represents a state 

outgoing to reactants, and for qi — >• +oo the wave function represents a 

state outgoing to products. 

Another pair of independent solutions of Eq. (j3.9ap . {qi\iljf'^^) and cor- 
responding to states "incoming from" reactants and products respectively, is obtained 
by requiring their momentum representations to be given by 



The corresponding position representations are obtained by way of Fourier transform, 



(gil-) = -=i= / e^'^-^^/'^^HpiDdp, , (3.16) 



of and (pi I "0/ defined in Eq. ()3.15p . namely 



cff Jo 



oo 







{qil^'"^) = ^=i= / " e^^iP^/^^ff (_p,)-i/2-^//?^o« dpi . (3.17b) 

The integrals in ()3.17p are not absolutely convergent, but can be defined as oscilla- 
tory integrals. The motivation for defining incoming states according to Eq. (|3.15p 
becomes clear from considering the stationary phase contributions to the integrals 
(j3.17p . These come from the pi satisfying the stationary phase condition 

-^(-/ln|pi|+(7ipi) = 0, (3.18) 
dpi 
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i.e., pi = I/qi, where pi > for and pi < for {qi\ip^j'^^ ) . This way one 

can associate the Lagrangian manifolds 



A 



(i,r) 



I 

(i,P) 



: Pl>0 



gi, — ) ■■ Pi<0}. 



(3.19) 



with the incoming states. These manifolds are also shown in Fig. One sees that for 
pi — )• +c« the wave function {qilipy ) represents a state incoming from reactants, 

and for pi — )• — oo the wave function represents a state incoming from 

products. Finally, calculating the Fourier integrals in Eq. ()3.17p [30j . one obtains 
expressions 



■ In h„ 



2 hcs 



«eff - ie ^ i^L \^p'f''^) \ , (3.20a) 



■ In tlr. 



2 Aleff 
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cff 1^/;} 



{o,p)> 



(3.20b) 



that relate the incoming states to the outgoing ones 
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Figure 7: The blue lines show the Lagrangian manifolds A^''^'*, A^''^\ A^°''^'*, and A^°''''* 

associated with the J-eigenstates I'?/'/ ''^■'), IV'/'^''), \'4^'f'^^)i and IV'/"'^^) respectively. The 
arrows indicate the classical Hamiltonian vector fields generated by / = (Pf — Q'\)/2. The 
red thick line corresponds to the dividing surface s(Qi, Pi) = Qi = 0. 

The equation for eigenstates of the action operator /, Eq. ()3.9ap . can also be 
solved in the Qi-representation [321 EDI- The two independent solutions are 

vr \2 hcsj 2+*A,fi \^ \l hes J 

(QilV'?'"^) = -(2neff)-i e^'^ r (i - i-/-) fe^'^./AQ,') , (3.21b) 

TT \2 n-eff/ ^ ''off V V ''-cfr / 
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where D,y is the parabohc cyhnder function of order u |51j . The outgoing states, 
and ars then obtained from Eq. ()3.20p . Fig. [7]shows the Lagrangian 

manifolds A^'^^^ A^''''', A^j''^\ and A^°''''' in the {Qi,Pi) coordinates. 

Finally, we note that the incoming and outgoing eigenstates of the / operator are 
orthogonal and satisfy the completeness relations |39l [50] 

\4'Wr^\ + \4''^){4''^\)dI = i, (3.22a) 
+ d/ = i , (3.22b) 

where 1 stands for the identity operator. 

Generally speaking, the qi- and Qi-representations of the eigenstates of the / 
operator are completely equivalent, and the choice of a representation is usually 
dictated by a particular problem. For example, in order to determine the total 
probability flux corresponding to the eigenstates through the dividing surface 

it is convenient to adopt the Qi-representation. Indeed, in the QP-basis the operator 
corresponding to the flux through the dividing surface s{Qi, Pi) = Qi is [30] 

F = -^[/,e(Qi)] = i(A<J(Qi) + 5(Qi)A) . (3.23) 



Thus, the expectation value of F with respect to the incoming states is given by 

d 

'dQi 



=0 



^ ^ (3.24) 



2Trhes 1 + exp(— 27r///ieff) 



The operators J^, defined by Eq. (j3.3p and corresponding to the action integrals 
of the bath modes, are the Hamiltonian operators of one-dimensional harmonic oscil- 
lators with unit frequency. Therefore, the wave function representing the eigenstates 
in Eq. (|3.9b[) are given by 

(I'^l^n,) = , f-^V^'^". exp (-4-^ , (3.25) 



V^2"'=(nfc!) V^^eff/ \VKsJ V 2/iefr 

where Hn denotes the nth order Hermite polynomial [51] . 

Here we note that the quantum normal form Hamiltonian can be interpreted as 
representing the activated complex. As it has been shown above the QNF Hamil- 
tonian has only continuous spectrum, and so there are no bound states associated 
with the activated complex. Physically this corresponds to the fact that the acti- 
vated complex has a finite life time. The latter is determined by the Gamov-Siegert 
resonances which we discuss in detail in Section [5j 
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4 The Cumulative Reaction Probability 
4.1 General Formulation 

The incoming and outgoing eigenstates of / defined in Sec. [3] are related to one 
another by Eq. (j3.20p . Therefore, each solution \ipj) of Eq. (j3.9ap can be written as 

a linear combination of 1^/°''^^'''*) or 1^/''^'^^^), 

li^j) = ap|4°'P)) + a.li^p'^) = /3p|vf 'P)) + /3r|Vf ''^> • (4.1) 
These representations are connected by the (one-dimensional) S-matrix, 



The entries of the S-matrix can be read of directly from Eq. (j3.20p yielding 



"eft 



S{I) = ^= r(--i— \ ■ (4.3) 

Using the relation r(l/2 + iy)r(l/2 — iy) = vr/ cosh(7ry) it is easy to see that 
S{I)*S{I) = 1 implying that the S-matrix is unitary. 

The transmission coefficient corresponding to the scattering matrix 5(1) reads 

T{I) = \Su{I)\' = = , (4.4) 

e ''off + e ''cff 1 + e ''cff 

and the reflection coefficient 

7^(/) = |5n(/)|' = / , = . (4.5) 

e ''eff + e off 1 + e off 

As required, one has T{I) +TZ{I) = 1. The characteristic action scale is given by the 
effective Planck's constant: T tends to 1 if I S> h^s and to if / <C —h^s- 

In accordance with Eqs. ()3.8p and ()3.9p . the incoming and outgoing scattering 
states of the d-dimensional reactive system at energy E are given by 

I^tCWp)) = l^'-'^/P)) ® IV'^) , (4.6a) 
l^rCo.r/p)) = |4°''-/P)) , (4.6b) 

where n = (77,2, • • • , n^) G INq^^ is a (d — l)-dimensional vector of scattering quantum 
numbers, and 

IV'n) = IV'nJ®...® IV'n,) (4.7) 

are the eigenstates of a (d— l)-dimensional harmonic oscillator corresponding to the 
bath modes. Here, / and n and not independent, but related by Eq. (|3.10p . For 
energies E close to the equilibrium point energy £^0 this equation can be resolved for 
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the reactive coordinate action, yielding a single valued function / = In{E). Then, 
the S-matrix connecting the scattering states |^(''''/p)) and |^'(°'''/p)) is block-diagonal 
and given by 

Sn,m{E)=6n,mS{In{E)), (4.8) 

where 5n,m is the multi-dimensional Kronecker symbol, and S{I) is given by Eq. ()4.3p . 

The transmission matrix T can be defined as a diagonal matrix with the diagonal 
elements equal to the modulus squared of the (1, 2)-components of the matrices S{I) 
in Eq. (gS]), i.e., 

Tn,m{E) = dn,^ |5i,2(/n(^))|' = ^n,: 

The cumulative reaction probability J\f(E) is then defined as (see, e.g., [3]) 

N{E) =Tt{T{E)] , (4.10) 
where Tr stands for the trace. Thus, using Eq. ()4.9p one obtains 

M{E) = Y,Tr.AE) = Y. 

n n 

The cumulative reaction probability N{E) is the quantum analogue of the clas- 
sical flux f{E) through the dividing surface or, more precisely, of the dimensionless 
quantity /(^)/(27r;ieff)'^"^ To see this one should consider M{E) in the semiclassical 
limit /loff —7- 0, where 

= G(/), (4.12) 

where is the Heaviside function. This means that the transmission coefficients 
Tn,n{E) in Eq. ()4.1ip are characteristic functions, i.e., in the semiclassical limit, 
Txi,n{E) is or 1 if In{E) is negative or positive respectively. This way the cumu- 
lative reaction probability can be considered to be a counting function. For a given 
energy E, it counts how many of the solutions /„ of the equations -K'QNF(-^ni ^efr(^2 + 
1/2), . . . , hcs{nd + 1/2)) = E with scattering quantum numbers n = (n2, . . . , n^) G 
Wq"^ are positive. In other words, M{E) can be considered to count the number of 
open "transmission channels" , where a transmission channel with quantum numbers 
n is open if the corresponding transmission coefficient T^^niE) is close to 1. 

Graphically J\f{E) can be interpreted as the number of grid points {hcs{n2 + 
1/2), . . . , heg{nd + 1/2)) in the space of action integrals (J2, . . . , Jd) £ [0, oo)'^""^ that 
are enclosed by the contour i^QNF(0, J2, ■ ■ ■ , Jd) = E, see Fig. [51 The number of grid 
points is approximately given by the volume in the space of (J2, . . . , Jd) G [0, (X))'^~^ 
enclosed by Kq^-f{0, J2, ■ ■ ■ , Jd) = E divided by ^eff*^"^- Using the fact that for 
^eff 0, -f^QNE becomes the function i^cNF which gives the classical energy as a 
function of the classical integrals (/, J2, . . . , Jd) we find that the volume in the space 
of (J2, . . . , Jd) enclosed by ircNF(0, J2, • • • , Jd) = E is given by the classical fiux f{E) 
divided by {2tt)'^~^, see Eq. (|2.14p . and the cumulative reaction probability M{E) is 



1 + exp 



2-K 



UE) 

^eff 



(4.9) 



1 + exp 
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UE) 



^eff 
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1 + exp - 2-Kl/h, 
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Figure 8: (a) Lines (/, hcs{n2 + 1/2), . . . , Hcsind + 1/2)), / £ M, rife e Fq, fc = 2, . . . , d, in the space 
(/, J2, • ■ • , Jd) G K, X [0, oo)'^~^ for d = 3 and their intersections with the surface i4rQNF(^, J2, J3) = E. (b) 
Grid points (/loff (?^2 + 1/2), • • ■ , hcsind + 1/2)) in the space (J2, . . . , Jd) for d — 3. The blue hue marks the 
contour iirQNF(0, J2, • ■ • , Jd) = E- In this plot only the scattering states for which the quantum numbers 
("•2, "-3) have the values (0,0), (0, 1), (1,0) or (1, 1) correspond to "open transmission channels", see text. 

thus approximately given by AAweyi(-E) = f{E)/{2TTh)'^-^ defined in Eq. (i2J5]l . This 
way we verify our statement in Sec. 12.41 that A/weyi (-E') gives the mean number of 
open transmission channels. In fact, as mentioned in Sec. 12.41 the classical flux f{E) 
can be considered to be the phase space volume enclosed by the energy contour 
of energy E of the invariant subsystem which has one degree of freedom less than 
the full scattering system and which as the so called activated complex is located 
between reactants and products. MweyliE) counts how many elementary quantum 
cells of volume (27r?iefr)'^~^ fit into this phase space volume and this way gives the 
Weyl approximation of the cumulative reaction probability N{E). 

It is important to note here that like the flux in the classical case the cumulative 
reaction probability is determined by local properties of the Hamilton operator H in 
the neighbourhood of the equilibrium point only. 

4.2 Nitrogen Exchange Reaction 

In this section we focus on the QNF theory as a tool for computation of the cumulative 
reaction probability (CRP) and the thermal reaction rate constant. To demonstrate 
the efficiency and capability of the QNF method and to compare it with other existing 
methods we follow Ref. [31j and address the CRP and the thermal rate in a collinear 
triatomic reaction, namely the nitrogen exchange reaction 

N + Ns^Na + N. (4.13) 

The collinear version of Eq. (|4.13|) corresponds to a two degree of freedom reactive 
system with the Hamiltonian operator of the form 

H = H{^i,^2,Pi,p2) = \ (pi+vi)+V{^i,^q2), (4.14) 
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where V{qi,q2) gives the Born-Oppenheimer potential energy surface, and iji and 
q2 are the Delves mass-scaled coordinates [52]. The effective Planck's constant is 
given by hgs = where n is the dimensionless reduced mass of the triatomic 

system (given in units of the electronic mass). For the case of nitrogen exchange 
reaction /x w 1.47 x 10^ yielding hcs ~ 8.2 x 10^'^. We further adopt the London- 
Eyring-Polanyi-Sato potential energy surface [S3| that possess a single saddle point 
governing the reaction from the asymptotic reactants and products states. 




E (eV) 

Figure 9: Cumulative reaction probability as a function of the total energy, M{E), for 
the coUinear nitrogen exchange reaction fl4.13p . The effective Planck's constant is hes ~ 
8.2 X 10~^. The J\f{E) curves obtained with the 4*^ and 6**^ order QNF are essentially 
indistinguishable for most of the energy range. The vertical dashed line shows the saddle 
point energy, Eq. 

Following the algorithmic procedure presented in Sec. [3] and Appendix |X] we 
compute the iVth order approximation 

^qnf = 4nf(^,'/2) (4.15) 

of the original Hamiltonian operator H. The CRP is then given by 

-1 

(4.16) 

with In^iE) solving 

Ki^M^, Ks{n2 + 1/2)) = E . (4.17) 



1 + exp —2-71 



In,{E) 



T-eflf 
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The resulting CPR-vs-energy curves computed for energies in the range between 1.3 
and 1.7 eV are shown in Fig. [9l The sohd color curves (blue, green, and red) cor- 
respond to different orders of the QNF approximation (2, 4, and 6 respectively). 
The N{E) curves obtained with the 4th and 6th order QNF are essentially indistin- 
guishable for most of the energy range; this fact signals the rapid convergence of the 
QNF expansion for the given value of the effective Planck's constant. The vertical 
dashed line shows the saddle point energy, £"0, of the London- Eyring-Polanyi-Sato 
potential energy surface. 

The circular data points in Fig. represent the "exact" valued of the CRP ob- 
tained through the full reactive quantum scattering calculation |54^ I55j . The lat- 
ter were performed by integrating the coupled multichannel Schrodinger equation 
in hyperspherical coordinates |54^ [55] from the strong interaction region to the 
asymptotic reactant and product configurations. The log-derivative matrix method 
of Manolopoulos and Gray [56] together with the six-step symplectic integrator of 
McLachlan and Atela [57] was used to integrate the radial Schrodinger equation. It 
is evident from Fig. [9] that the quantitative agreement of the exact and QNF values 
of N{E) well extends up to energies of 1.5 eV. d 

We finally note that the QNF calculation of the CRP requires significantly less 
computational time than the corresponding full quantum reactive scattering calcu- 
lation. For example, the 6th order QNF computation of the nitrogen-exchange CRP 
curve in Fig. [Uj took about 10 minutes on a 2.6 GHz processor, 2 GB RAM com- 
puter, while the corresponding full quantum reactive scattering computation took 
more than 12 hours on the same machine. The QNF approach becomes even more 
advantageous for treating chemical systems of atoms heavier than nitrogen: the ex- 
pense of the full quantum computations rapidly grows with the number of asymptotic 
channels (and, therefore, with mass) [58], while the QNF expansion only becomes 
more rapidly convergent making the corresponding analysis computationally cheaper. 

We now discuss the straightforward extension of the QNF approach to computa- 
tion of the thermal reaction rate constant, fc(T), defined as [59] 



where T stands for the absolute temperature, Qt{T) is the partition function of the 
reactant, and /cb is the Boltzmann constant. To this end, we use the CRP data for 
the nitrogen exchange reaction, Eq. ()4.13p . presented in Fig. [9l The data allow us 
to compute k{T) in the range of temperatures between 600 and 700 K, for which 
the integrand in the right hand side of Eq. (|4.18p is well localized to the energy 
interval between 1.3 and 1.7 eV. Figure fTOl provides a comparison of k{T) calculated 
from N{E) obtained using the QNF of orders 2, 4, 6 (color solid lines), and the 
"exact" , reactive quantum scattering technique (black line with circles) . The inset in 



^It is interesting to note that the QNF curves in Fig. IH] seem to converge to a result that slightly differs 
from the exact CRP curve at high energies: a small discrepancy starts to show up at energies around 
1.5 eV and becomes more pronounced at energies above 1.65 eV. This discrepancy may be attributed to 
the growing importance of corner-cutting tunneling at high energies. We discuss the phenomenon of the 
corner-cutting tunneling in relation to the normal form theory in Section [6. II 
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Figure 10: The thermal rate constant (multiphed by the reactant partition function) as 
a function of temperature for the colhnear nitrogen exchange reaction fl4.13p . The curves 
obtained with the 4**^ and 6**^ order QNF are essentially indistinguishable for most of the 
energy range. The inset shows the ratio, ^qne/^qs, of the thermal rate constant computed 
using the QNF of orders 2, 4, 6, to the one obtained from the quantum scattering data. 



the figure shows the ratio, A;qnf/^qs> of the thermal rate computed using the QNF 
of orders 2, 4, 6, to the one obtained from the quantum scattering data. One can 
clearly see that the thermal rate constant computed with the QNF method rapidly 
approaches its exact value as the approximation order N is increased. 



4.3 A 3 DoF Reactive System 

In order to illustrate the utility of the QNF technique for computation of reaction 
rates in higher dimensional systems we consider a 3 DoF model system consisting of 
an Eckart barrier in the (physical) gi-direction that is coupled to Morse oscillators 
in the g2-direction and in the gs-direction. The Hamiltonian operator is 

H = \{M+M+Pl) + ^e(Ii) + Vm;2{I2) + Vud^s) + eHc ■ (4.19) 



Here, 



( ) = ^ exp((g + go)/a) ^ exp((, + ,o)/a) 

l + exp(((? + go)/a) [1 + exp((Q + (/o)/a)]2 ' 
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with positive real parameters a, A, and B, and 

, B + A 

is the Eckart potential. For B > A > it possesses a maximum at g = 0; the 
potential value at the maximum is Ve(0) = {A + Bf/{AB). The Morse potentials, 
given by 

yM;k{q) = De-k [exp (-2aM;fc q) -2 exp {-au-k q)] , (4.22) 

are parametrized by -De;fc and aM;fc, for A; = 2,3 respectively. For He we choose the 
mutual kinetic coupling 

He =PlP2 + P2 P3 +P3P1- (4.23) 

The strength of the coupling is controlled by the parameter e in Eq. (j4.19|) . The 
vector field generated by the corresponding classical Hamilton function has an equi- 
librium point at {qi,q2,q3,Pi,P2,P3) = 0. For |e| sufficiently small (for given values 
of parameters of the Eckart and Morse potentials), the equilibrium point is of saddle- 
centre-centre stability type. 

In the uncoupled case, e = 0, it is easy to calculate the exact CRP, Mexa.ct{E), of 
the Eckart-Morse-Morse system analytically. Indeed, in accordance with Eq. (|4.1ip . 

oo oo 

■f^exactiE) = ^ ^ ^exact (E — EM;2,n2 — -E'M;3,n3) , (4.24) 
"■2=0 n3=0 

where Texact denotes the exact transmission coefficient for the Eckart barrier, and is 
given by [60] 

cosh(27r(a - /3)) + cosh(27r(5) 
cosh(27r(Q; + pj) + cosh(27ro) 

with 



1^1 I E- A 1 I B-C h,^ 

"=2VC' ^=2V^' ^=2V^' ^ = ^- ^'-^'^ 

In Eq. (j4.24p . Ey[.k,nk (with /c = 2,3) are the energy levels of the one-dimensional 
Morse oscillators. 



1 2 2 1 xf^De-k \ 

EM;k,nk = --^O-M-kf^cS \ nk + , Tl-fc = 0, 1, 2, . . . . (4.27) 

2 ' y Z aM;k fteff J 

The graph of Afexact{E), obtained from Eqs. ()4. 24114. 27ll . is shown by the black oscil- 
latory curve in the top panel of Fig. [TTJ Here, the parameters for the Eckart potential 
are a = 1, A = 0.5, and B = 5, the parameters for the Morse potential are De-2 = 1, 
De-3 = 1.5, and aM-2 = clm;3 = 1, and the effective Planck's constant h^s = 0.1. 

The quantum normal form computation of the CRP for the Eckart-Morse-Morse 
system gives 



^f<iME)=Y: E [l + expf-2vr%^) 



(4.28) 
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Figure 11: The top panel shows the cumulative reaction probabilities J\fexa.ct{E) (black 
oscillatory curve) and A/weyi(-E') (red smooth curve) for the Eckart- Morse- Morse reactive 
system with the Hamiltonian given by Eq. f l4.19p with e = 0. It also shows the quantum 
numbers {n2,n3) of the Morse oscillators that contribute to the quantization steps. The 
bottom panel shows the resonances in the complex energy plane marked by circles for the 
uncoupled case e = and by crosses for the strongly coupled case e = 0.3. The parameters 
for the Eckart potential are a = 1, A = 0.5, and B = 5. The parameters for the Morse 
potential are De;2 = 1, -De;3 = 1.5, and aM;2 = ciM;3 = 1- Also, hes = 0.1. 



with In^^nsiE) solving 



^qnf(^. ^cfr(n2 + 1/2), ?ieff(n3 + 1/2)) = E . (4.29) 

The high quality of the quantum normal form approximation of the cumulative reac- 
tion probability is illustrated in Fig. W2\ which shows IA/qnf — AAexact I as a function 
of the energy E for different order N of the quantum normal form. It is evident that 
the CRP computed with the QNF method rapidly approaches its exact value as the 
approximation order is increased. 
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Figure 12: Errors |A/qnf — A/'exactI for the cumulative reaction probability as a function of 
energy E of the Eckart-Morse-Morse system for different order N of the quantum normal 
form. 

5 G amov- Sieger t Resonances 

As we have discussed in Section [2] in classical mechanics the transition state is rep- 
resented by a lower dimensional invariant subsystem, the NHIM. In the quantum 
world, due to Heisenberg's Uncertainty Principle, we cannot localise quantum states 
entirely on the NHIM, so there can not be any invariant quantum subsystem rep- 
resenting the transition states. Instead we expect a finite lifetime for the transition 
state. The lifetime of the transition state is determined by the Gamov-Siegert reso- 
nances, whose importance in the theory of reaction rates has been emphasized in the 
literature |^ i62j . 

Intuitively, a resonance is a complex eigenvalue of a Hamiltonian, 

H\i;res) = E\i^res) (5.1) 

where the imaginary part has to satisfy ImE < 0. The time evolution of a resonance 
state is 

and so we see that a resonance state decays in time with lifetime given by 

~ \lmE\ 

Notice that since H is self-adjoint, a resonant state with complex energy cannot be 
square integrable. 
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In order to see how the resonances can be computed from the QNF Hamiltonian, 
consider first the example of a one-dimensional system described by the Hamilton 
operator 

i = ^iqp + pq)- (5.2) 

We easily see that |n) defined by {q\n) = is a resonance eigenstate 

i\n) = -ihefiin + l/2)\n) . (5.3) 

with complex eigenvalue 

En = -ihea{n + 1/2) , n = 0,...,N. (5.4) 

These are the resonances of the operator /. Notice that compared to the scattering 
states |/) which depended on a continuous parameter /, the resonance states are 
quantized and depend on a discrete parameter n € Mq . Using this simple example we 
can now directly determine the resonances and the corresponding resonance states for 
a Hamilton operator in quantum normal form Hq^f = Kq^p{I, J2, . . . , Jd)- Let ipn,, 
denote the n)^ harmonic oscillator eigenfunction (see (|3.25p ). For n = (ni, . . . , n^) G 
IN^, set 

\ni,n2, ■■■ ,nd) = \ni) (g) {ipn^) iX" • • • l^n^) , (5-5) 

then we have 



-f^QNF|"-l,n2, • • • ,nd) 

= Kqj^f{ - ihesirii + l/2),;ieff(n2 + 1/2), . . . ,Ks{nd + 1/2)) |ni, n2, • • • , 

(5.6) 

and so the resonances of -ffqNF are given by 

En^,n2,...,ni = i^QNF ( " iKs{ni + 1/2), Ks{n2 + 1/2), . . . , ^efr("-d + 1/2)) . (5.7) 

The bottom panel of Fig. [11] shows the resonances in the complex energy plane that 
have been computed using Eq. (|5.7|) for the case of the coupled Eckart-Morse- Morse 
reactive system discussed in Section 14.31 

The resonances computed from the QNF describe the lifetime of the activated 
complex. To see this in more detail consider a state localized at time t = on the 
NHIM, i.e., the dependence on the local normal form reaction coordinate qi is of the 
form ^ 

(gi|yo)= , .\i/4 e~^^ • (5.8) 

which is a minimal uncertainty state. In the bath coordinates we take the state to 
be given by harmonic oscillator eigenfunctions V'nfe, so that 

1^) = IVJO)® IV'na)®---® IV'n,) (5.9) 

for some fixed quantum numbers n2,---,nd G INq. Then, expanding |^') into the 
basis of Eq. (|5.5|) we obtain the time evolved state 

= (^^S^E ^^(ki^''-*"- -"'■•|2n,n,.. . .,„,). (5.10) 
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A suitable quantity for measuring the lifetime of a state like (|5.9p is the decay of 
the autocorrelation function 

m^m'- (5.11) 

For the Hamiltonian in quantum normal form the overlap reads 

(^l^(t)) = ('iy^'^Il!i±i/^(_i)"e-*^^2„,„2,...,„/Meff , (5,12) 

n=0 

The leading term in this sum for t — )• oo is given by the smallest resonance with 
n = 0. Hence, 

|(^|^(t))|2 ~ 2e2i^-^«'"2>--"d*/'"'-ff , (5.13) 

and this determines the maximal lifetime of a quantum state of the activated complex, 
i.e. a state initially localised near the classically invariant subsystem given by the 
NHIM. 

For small h^g the quantum normal form is dominated by its quadratic part and 
that gives 

2Im^o,n2,...,nd/^eff ~ -A, (5.14) 

and therefore for small ?ieff 

|(^|^(t))|2 ~2e-*^ . (5.15) 

Thus, the quantum lifetime of the activated complex is in leading order for small h^s 
given by the reciprocal value of the classical Lyapunov exponent associated with the 
saddle equilibrium point. 



6 Further Challenges 
6.1 Corner-Cutting Tunneling 

The approach based on normal forms is designed to give an accurate description of the 
dynamics near the saddle, or near the transition state. If the reaction is influenced by 
processes which avoid the neighborhood of the saddle then the normal form is unlikely 
to capture them. One such process is corner cutting tunneling. We want to point 
out that the results derived from the quantum normal form, like the computation of 
the S-matrix and cumulative reaction probability, contain tunneling contributions, 
but these are related to tunneling paths which are close to the saddle. One strength 
of our approach is that it describes uniformly the transition from energies below the 
saddle energy, where only tunneling is possible, to energies above the saddle energy. 

But there are situations where other tunneling contributions have a large effect. If 
the reaction path (e.g., the minimum energy path) in configuration space is strongly 
curved in the regions where it connects reactants and products, then the reactant 
valley and the product valley can have a small angle relative to each other. In 
this situation there can be quite short tunneling paths which connect reactants and 
products and which don't go near the saddle region. So these tunneling paths "cut 
the corner". The contribution from such trajectories play an important role, in 
particular for hydrogen exchange reactions at room temperature and lower [63]. 
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In these situations the normal form approach has to be combined with multidi- 
mensional semiclassical tunneling calculations. 

6.2 State-to-State Reaction Rates 

It is important to emphasize again that, so far, our approach to quantum reaction 
dynamics has been local, i.e., it is derived completely from the properties of the 
quantum normal form that is valid in the neighborhood of the saddle-centre-- • • - 
centre equilibrium point. The property of the resulting S-matrix in ()4.8p being block- 
diagonal reflects the fact that the quantum normal form is integrable in the sense 
that the basis of scattering states can be chosen in the product form ()4.6p . In a 
different basis the matrix will lose this feature, and phenomena like mode mixing 
are related to how other incoming and outgoing scattering states are related to this 
special basis. It is natural to embed the study of this phenomenon in a study of the 
global dynamics which we will describe in this section. The global formalism is in 
particular required in order to compute general state-to-state reaction rates. 

Let us start by describing the scattering or reaction process in classical mechanics 
by using Poincare sections. Recall that a Poincare section at energy E is given by a 
smooth hypersurface S(-E') of the energy surface with energy E which is transversal to 
the flow is allowed to have several components). If we have two such Poincare 

sections Tii{E) and Ti2{E) such that all the flow lines intersecting T,i(E) intersect at 
a later time T,2{E), too, then moving along the flow from T,i{E) to Ti2{E) defines a 
Poincare map 

p(2,i)(^):Si(E)^S2(^) . (6.1) 

Such Poincare maps can be composed. If T,3{E) is another Poincare section which lies 
behind S2(-B) in the sense that the flow lines that intersect S2(-E) also intersect ^^{E) 
at a later time, and if P^^''^\E) : T,2{E) — )• T,-i{E) is the corresponding Poincare map, 
then the Poincare map 

p(3,i)(i^):Si(i^)^S3(i^) (6.2) 

is given by 

p(3,l) = p{3,2) ^ p(2,l) _ 3) 

Using this construction we can describe transport through phase space regions 
by a sequence of maps. Given some Poincare section Sinitiai(i^') located in the area 
of initial points in the reactants region where we prepare the system and a Poincare 
section Ilfinai(-£') in the products region where we measure the outcome, a succession 
of Poincare maps 

Sinitial(^) ^ Si(^) ^ S2(^) ^ > Sfinal(^) (6.4) 

tells us how the initial points are transported through the systemj^l 

The advantage of subdividing the flow into a sequence of maps lies in the fact 
that different regions in phase space might need different techniques to compute the 

^ We here ignore the difficulties involved in constructing global Poincare sections (see, e.g., [64]); we 
assume that the sequence of Poincare sections (|6.4p is intersected transversally by the trajectories with 
initial points from a suitable open subset in the reactants region. 
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flow. In our case of interest Poincare sections can be constructed to the products and 
reactants side of a saddle-centre-- • • -centre equilibrium point. The dynamics 'across' 
this equilibrium point can then be described by the normal form while the dynamics 
between neighbourhoods of different saddle points can be obtained from integrating 
the original equations of motions [651 |66l HZ] • Moreover, the phase space structures 
obtained from the local normal form can be "globalized" following the discussion in 
Section 12.61 

A similar procedure can be developed in the quantum case. The Poincare maps 

(6.5) 

are symplectic maps, and as such can be quantised using the theory of Fourier integral 
operators. The quantisations will be unitary operators which we interpret as local 
S-matrices, 

S^^'^He) : Ll^j,^ ^ Ll^^E^ , (6.6) 

where L'^^^^-^ is a Hilbert space obtained by geometric quantisation of ^{E), see, e.g., 
[67j . This is similar to the quantisation developed in [68] . As in classical dynamics 
we can compose these matrices to obtain a global S-matrix 

^(final, initial) g{{inal,n) ^j^-^ g{n,n~~l) ^jjj-j _ _ _ ^(1, initial) ^^■^ 

which tells us how initial states in Li. are transformed into final states in 

i^initiaU-t-j 

^Sfl i{E)' "^^^ reasons for introducing this splitting of the S-matrix are the same as in 
the classical case. We can employ different techniques for computing the S-matrices 
according to different local properties of the system. Near equilibrium points the 
dynamics can be described by the quantum normal form we developed in this paper. 
Notice that the neighbourhoods of the saddle-centre-- • • -centre equilibriuml points 
are the regions where we expect quantum effects to be of most importance due 
to partial reflection at and tunnelling through the barriers associated with saddle 
points. The quantum transport between neighbourhoods of different equilibrium 
points can be described by a standard van Vleck type formalisms, using, e.g, initial 
value representations (IVRs) which are very common in theoretical chemistry (see, 
e.g., [11[69] for references). 

6.3 Flux-Flux Autocorrelation Function Formalism 

In this section we discuss the flux-flux autocorrelation function approach to compu- 
tation of reaction rates in both classical and quantum theories. 

Classical formulation: The directional flux through the dividing surface that 
determines the classical reaction rate, see Sec. 12. 4| can be written as |8l [9l 0] 



f{E)= [ 5{E-H{z))F{z)P,{z)dz. (6.J 

yiR2d 



The dynamical meaning of Eq. (|6.8p is as follows. First, z = {qi,Pi, ■ ■ ■ ,qd,Pd) is 
a point in the system's phase space, and the function 6{E — H{z)) restricts the 



39 



integration to the energy surface of energy E under consideration. The flux factor F 
is given by 



(6.9) 



where zt is the Hamiltonian flow satisfying the initial condition zt=o = z. The dividing 
surface is defined as the zero level set of the function s, i.e., {z G IR-^'^ : s{z) = 0}. It 
is assumed that this surface divides the phase space into two components: the region 
given by s{z) < corresponds to reactants, while the region given by s{z) > to 
products. Finally, the projection function is defined as 

P,{z) = lim e{s{zt)) , (6.10) 

and equals unity (zero) if a trajectory starting at z ends up in the product (reactant) 
region for infinitely long times. 

In order to explain the relation of Eq. (|6.8p to Eq. ()2.14p we first rewrite the 
projection function as 

(•oo 



/■oo J 

PAz) = e{s{ztj) + -e{s{zt))dt, (6.11) 



where we choose to = ~e with e — )• 0^. At this stage it is important to note that since 
the flux factor F{z) is proportional to 6[s{z)) the integral in Eq. (jG.Sp is effectively 
restricted to the dividing surfacJ^. But, using such properties of the dividing surface 
that it is (i) recrossing-free and (ii) transverse to the Hamiltonian flow, we have 

lim Q(s(z^) = 1 - PJz) (6.12) 

for z on the dividing surface. Then, substituting Eqs. (|6.9p and (|6.12p into Eq. (|6.1ip . 
and solving for the projection function Pr we obtain 

P,(z) = -+ hm - / F{zt)dt (6.13) 

for z on the dividing surface. Finally, substituting Eq. ()6.13p into Eq. ()6.8p . and 
taking into account the fact that the total flux through the dividing surface is zero, 
i.e., 

/ 5{E - H{z))F{z)Az = {), (6.14) 
we arrive at the expression for the directional flux, 

f[E)= \ r CF{t)dt, (6.15) 



^"More precisely, since in addition to the flux factor F one also has to take into account the function 
d{E — H) the integral in Eq. (|6.8p is a {2d — 2)-dimensional integral over the intersection of the dividing 
surface with the energy surface of energy E. 



40 



as a time integral of the flux-flux autocorrelation function 



CF{t)= ! 5{E - H{z))F{z)F{zt)dz . (6.16) 

We now explicitly calculate the time dependence of the function Cp using the 
classical normal form theory. To this end we first express the flux factor F in the 
normal form coordinates (Qi , Pi , 52 j P2 , • • • , QdiPd) , see Sec. 12. H in which the dividing 
surface is given by s{z) = Qi: 



F(z) = ^e(Qi(t)) 



t=o 



(6.17) 



t=o 



where 



d 

A = -qjHcnf{I, J2, ■ ■ ■ ,Jd) 



(6.18) 



and -f^cNF is the normal form Hamiltonian as a function of the action integrals 
/ = (Pf - Ql)/2 and = {pi + ql) /2, k = 2, . . . , d. Equation (pT^ then reads 



CF{t)= [ 6iE-HcNF)6iQi)d{Qiit))A^PiPi{t)dz. (6.19) 

The product of the 5- functions of Qi and the corresponding time evolved coordinate 
Qi{t) signals that only the infinitesimally short time scales t — )• give a non- vanishing 
contribution to the integral. Linearization the time evolution of zt around t = yields 

Cpit) = 6{t) [ 6{E - i/cNF)5(Qi)A|Pi| dz 

= 25{t)i27r)^-' [ 5{E - Fcnf)^^^ d/dJs . . . d 

= 25{t){2T:f~^ [ dJ2...dJrf. (6.20) 

The last integral in Eq. ()6.20p is nothing but the volume V{E) in the action space 
(J2,...,Jd) enclosed by the contour Hcnf{0, ■ ■ , Jd) = E, and illustrated in 
Fig. [5l that defines the directional flux through the dividing surface, see Eq. (j2.14p . 
Thus, we obtain 

Cpit) = 26{t){2TTy-^V{E) = 26{t)f{E) , (6.21) 

which is in agreement with Eq. ()6.15p . Finally, we note that in view of Eq. (|6.2ip 
one can rewrite Eq. (|6.15p as 



1 r+co 

f{E) = - j ^ CF{t)dt. 



(6.22) 
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Quantum formulation: In the quantum version of the flux-flux autocorrelation 
function approach Eq. ()6.8p for the directional flux is replaced by a corresponding 
equation for the cumulative reaction probability, 

M{E) = iTTh^s T> {S{E - H)FPr^ . (6.23) 



Here, the flux factor operator F is given by 

i 



t=0 ^eff 



[H,eis)], (6.24) 



where G(s) is a quantization of the composition of the Heaviside function with a 
function s defining the dividing surface. The projection operator Pj- is defined as 

Pr = lim e'^^/^-«e{s)e-'^^/^-^ . (6.25) 

Then, as shown in Ref. [4j, Eq. ()6.23p can be written in the close analogy with its 
classical counterpart Eq. (j6.15|) : 

M{E) = (2^/ieff) - / Cpit) dt , (6.26) 

where 

Cpit) = TV |(5(^ - F)Fe*^*/''-ffFe-*-f^*/'"'°ff } (6.27) 
is the quantum flux-flux autocorrelation function. 

As we have shown above the classical flux- flux autocorrelation function Cp can 
be explicitly analyzed within the framework of the classical normal form theory. The 
following natural question arises: can one obtain the time dependence of the quantum 
flux- flux autocorrelation function Cp using the methods of the quantum normal form 
theory? Of course, the QNF technique provides one with an approximation of the 
original Hamiltonian operator. This approximation is only accurate in the vicinity 
of the equilibrium saddle point in phase space, so one should not expect a perfect 
agreement between the QNF flux-flux autocorrelation function and the exact one to 
hold up to infinitely long times. Instead, the QNF theory will provide an approxima- 
tion of the exact flux-flux autocorrelation function in a certain time interval whose 
length will depend on the effective Planck's constant among other parameters. 



6.4 Convergence of Quantum Normal Form 

Both the classical and quantum normal forms have the form of power series where 
each successive term is constructed via an iterative technique. Consequently, in 
general, we will only be able to compute a finite number of terms of the normal form. 
Therefore the obvious question that arises is how many terms of the normal form 
are required in order that the quantities derived from the normal form are accurate? 
A discussion, as well as some references with specific examples, was given in Section 
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12.51 for the classical normal form. In this section we are concerned with the behavior 
of the quantum normal form, for which there has been essentially no work from this 
point of view. 

First, we begin with a brief background discussion. Since the normal forms, 
both classical and quantum, are given as series representations of the Hamiltonian 
function and Hamilton operator, respectively, a natural question to ask is do these 
series converge? In the classical setting there has been previous work on this question. 
However, practically speaking, from the point of view of using the results of the 
theory, there are three related questions: 1) convergence of the series representing 
the normal form Hamiltonian, 2) convergence of the transformation from the original 
coordinate to the normal form coordinates, and 3) convergence of any integrals of 
the motion that arise from the normalization procedure. A recent discussion of 
these issues, as well as a discussion of earlier results for the classical Hamiltonian 
normal form setting can be found in [TO]. Briefly, the situation with respect to 
convergence for Hamiltonian systems with three or more degrees-of-freedom is not 
optimistic. Generically, divergence is the expected behavior in normal form theory 
for classical Hamiltonian systems. However, the situation is not so pessimistic as 
one my initially believe. There are many examples of divergent series which still 
yield useful information in an asymptotic sense. Resummation techniques and Fade 
approximation techniques can be used to find an optimal number of terms that yield 
a desired accuracy. For classical Hamiltonian normal forms these issues have been 
examined in [Til [72l [73l [71]. These issues have yet to be explored in the quantum 
normal form setting. 

It is useful to note that the case of two degree-of-freedom Hamiltonian systems 
is special. In this case a classical result of ^75j (see also [76j) gives convergence 
results for the classical Hamiltonian normal form in the neighborhood of a saddle- 
center equilibrium point. Recently, the first results on convergence of the quantum 
normal form have appeared. In [77j convergence results for a one and a half degree- 
of-freedom system (i.e. time-periodically forced one degree-of-freedom Hamiltonian 
system) have been given. It is not unreasonable that these results can be extended to 
the quantum normal form in the neighborhood of a saddle-center equilibrium point 
of a two degree-of-freedom system. 

In the following we provide a qualitative discussion of the convergence of the QNF 
based on our calculations performed for the triatomic collinear reactions. In this 
situation the QNF approximates the Hamiltonian of the reacting system in a phase- 
space neighborhood of the saddle-center equilibrium point. Thus, for instance, in 
computing the CRP one only expects this approximation to render reliable results in a 
certain energy range around the saddle point energy Eq of the potential energy surface 
under consideration. The energy difference {E — Eq) may therefore be considered as 
one small parameter in the QNF expansion. The role of the other small parameter 
is played by the effective Planck's constant, fieff- It is the convergence of the QNF 
with respect to this second small parameter that we focus on in this section. 

We proceed by considering the right hand side of Eq. (j3.ip . i.e., the QNF, at / = 0, 
corresponding to no "energy" in the reaction coordinate, and n2 = 0, giving the 
zero-point 'vibrational energy' of the transverse degree of freedom. Then, Eq. (jS.ip 
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becomes 



[N/2\ 

E = Eo+ Cn^eff" . (6.28) 

ra=l 

For the case of the cohinear hydrogen exchange reaction, H + H2 — )• H2 + H, on 
the Porter-Karplus potential energy surface the first five expansion coefficients were 
obtained in Ref. [31]: ci = 0.161982, C2 = 1.193254, C3 = 14.90023, C4 = 378.7950, 
and C5 = 1227.035. As — )• 00 the radius of convergence /leff^"'' of the sum in 
Eq. (j6.28p is given by 

neff(°) = Mm ^ . (6.29) 

Here, we make a crude estimate of ^eff^*^^ by only considering the first five expansion 
coefficients in Eq. (|6.29p . i.e., c„ with n = 1, . . . , 5; then, the radius of convergence 
is given by ^eff*-^'' ~ 0.04. 
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Figure 13: Cumulative reaction probability as a function of the total energy, M{E), for the 
collinear hydrogen exchange reaction. Figure (a) corresponds to isotope characterized 
by ^eff = 3.07 X 10~^, (b) corresponds to isotope characterized by ^efr = 1-77 x 10~^, 
and (c) corresponds to a hypothetical "isotope" characterized by /Zeg = 6.9 x 10^^. 
The vertical dashed line shows the saddle point energy, Eq. 

The estimated value of h^^^'^^ sheds light on the seeming inefficiency of the QNF 
theory for CRP computations in light atom reactions. Indeed, the exchange 
reaction, see Fig. [Oh . is characterized by h^Q = 3.07 x lO^'^. This value being close 
to ^eff^^^ signals that the corresponding QNF expansion converges very slowly, if at 
all, and, possibly, terms of orders far beyond = 10 are needed for a reliable CRP 
prediction in Fig. [T3b . 

In the case of the exchange reaction the effective Planck's constant is h^s = 
1.77 X 10~^ and is thus smaller than h^^^'^^ This fact is in agreement with the 
apparent speed-up of the convergence of the CRP values, see Fig. [T3b . in comparison 
with the case. Finally, the convergence is very fast and pronounced for the case of 
the heavy (hypothetical) '^'^H atoms, see Fig. [T3b . for which = 6.9 x 10^'^ which 
is much smaller that the estimated convergence radius. 
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Clearly, there is a great deal of scope for further studies of convergence aspects 
of the quantum normal form as well as related optimal truncation and resummation 
techniques. 

7 Conclusions 

We have reviewed a new dynamical systems, phase space approach to quantum reac- 
tive scattering. The quantum theory arises from an underlying classical theory that 
reveals the geometrical structures in phase space that govern trajectories evolution 
from reactants to products in the reaction region. The classical theory is valid for 
any number of degrees-of-freedom, and the new concept which has led to this is the 
introduction of the notion of a normally hyperbolic invariant manifold, or NHIM, 
into classical reaction dynamics. The NHIM is the key phase space structure that 
leads to the construction of a surface dividing reactants from products having the no 
(local) re-crossing and minimal flux properties. Moreover, the stable and unstable 
manifolds of the NHIMs have one less dimension than the energy surface (i.e. they 
are "co-dimension one" ) and they bound regions of the energy surface (being invari- 
ant manifolds, or "impenetrable barriers") which contain all reacting trajectories. 
These phase space structures area realized in specific systems through the use of 
the classical Poincare-Birkhoff normal form. Quantization is carried out through an 
analogous quantum normal form theory and the Weyl quantization procedure. The 
related structure of the classical and quantum normal form theories makes the quan- 
tum manifestation of the classical phase space structures transparent. The phase 
space structures in the classical case define the "landscape" in the energy surface 
which constrain the location and evolution of reacting trajectories, and this renders 
the need for the calculation of trajectories to evaluate, e.g. classical flux across the 
dividing surface, unnecessary. This characteristic is inherited in the quantum setting 
in the sense that quantum expressions governing reaction are expressed in terms of 
quantities that can be computed from the quantum normal form without the need 
to compute classical trajectories. Thus this phase space approach to quantum reac- 
tive scattering provides a completely new approach to the computation of quantities 
describing quantum reactive scattering which may prove to be fruitful in the study 
of "large" quantum systems. 

Appendix 

A The Normal Form Algorithm: Classical and 
Quantum 

In Sec. El we mentioned that the crucial point to realize the phase space structures 
governing reaction dynamics is to choose a suitable set of phase space coordinates 
in terms of which the Hamilton function assumes a very simple form. Similarly, the 
description of the activated complex in Sec. [3] leading to the explicit expressions for 
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the cumulative reaction probabilities and the Gamov-Siegert resonances as explained 
in Sections 14.11 and \5\ respectively, relied on expressing the Hamilton operator in 
a suitable basis. Both the canonical (or equivalently symplectic) transformation of 
the phase space coordinates and the corresponding transformation of the Hamilton 
function in the classical case, and the unitary transformation of the basis set together 
with the corresponding transformation of the Hamilton operator in the quantum 
case can be constructed in an algorithmic fashion. The algorithms are based on a 
(classical) normal form (CNF) in the classical case and on an analogous quantum 
normal form (QNF) in the quantum mechanical case. In the following we give a 
brief description of the algorithms to compute these normal forms. We describe 
the algorithms in such a way that the similarities between both algorithms become 
apparent. For rigorous mathematical statements, proofs and further details we refer 
to Ref. [291 ED]. 

A.l The classical normal form algorithm 

We begin by considering a Hamiltonian system with phase space x with phase 
space coordinates q = {qi,q2, ■ ■ ■ , qd) and p = (pi,p2) • • • iPd) which we group in the 
2d-dimensional vector z = (q, p), and a Hamilton function (q, p). For convenience, 
we can use atomic units. This imposes no restriction. However, it has the advan- 
tage that the phase space coordinates become dimensionless. We will now assume 
that Hamilton's equations associated with -ff (q, p) have a (single) equilibrium point, 
zq = (qO)Po)) of saddle-center-. . .-center stability type. By this we mean that the 
matrix associated with the linearization of Hamilton's equations about this equilib- 
rium point has two real eigenvalues, itA, of equal magnitude and opposite sign, and 
d—1 purely imaginary complex conjugate pairs of eigenvalues ±iwfc, k = 2, . . . ,d. If 
the Hamilton function is of the form kinetic energy plus potential energy then this 
type of equilibrium point of Hamilton's equations correspond to an index one saddle 
point of the potential energy. 

The CNF procedure consists of a sequence of canonical (or equivalently sym- 
plectic) transformations changing the phase space coordinates z = (q, p) in a phase 
space neighborhood of the equilibrium point zq. This is accomplished order by order 
according to the sequence 

i ^ z(O) ^ ^(1) ^ ^(2) ^ ^(3) ^ . . . ^ ^(N) ^ (A.l) 

where z^") is obtained from z^"^^^^ by means of a symplectic transformation 

^ = </,vK„z^""^^ . (A.2) 
generated by a homogenous polynomial Wn{z) of order n, i.e. 

Wn G W3 := spanjg;^^ . . . q'^^p^^' ...p^/ : \a\ + |/3| = n} . (A.3) 

More precisely, the 0m/„ hi ()A.2p denote the time-one maps of the flows generated by 
the Hamiltonian vector fields corresponding to the polynomials Wn (see [30] for the 
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details). The maximum order N in (jA.ip is the desired order of accuracy at which 
the expansion will be terminated and truncated. 

Expressing the Hamilton function H in the coordinates z^"), n = 1,. . . ,N, we 
get a sequence of Hamilton functions H^"''' , 

H ^ H^o) ^ ^ ^ i/(3) ^ . . . ^ i^(^) , (A.4) 

where forn = 1, . . . , A^, and i7(")(z('^)) = i/("-i)(z("-i)) = ii'("-i)((/.^:^^z(")), i.e. 

i^W oc/,^^. (A.5) 

To avoid a proliferation of notation we will in the following neglect the superscripts 
(n) for the phase space coordinates. 

In the first transformation in (jA.ip we shift the equilibrium point zq to the origin, 
i.e. z I—)- (j)]y-^{z) := z — zq. This gives 

i/W(z) = f(°)(z + zo). (A.6) 

The next steps of the CNF algorithm rely on the power series expansions of , 

oo 

H^^\z)=Eo + ^H("\z), (A.7) 

s=2 

(n) 

where the Hg are homogenous polynomials in W^j given by 

ai!...a.!/3i!.../3,!^i ' ' ' •••^^ ^^'^^ 

with ak,(3k G Nq, |a| = Y^f^ak, = X]fc/3fc. For n = 1, the coefficients in (|A.8|1 are 
given by the Taylor expansion of iJ(^) about the origm, 



- fr ^^F«(z) 



(A.9) 



z=0 



For n > 3, the coefficients in (jA.Sp are obtained recursively. For n = 2, i.e. the 
second step in the sequence of transformations (jA.ip . the coefficients in (jA.Sp are 
determined by a linear transformation of the phase space coordinates according to 

z (pWii^) ■■= Mz. (A.IO) 

Here, M is a symplectic 2(i x 2(i matrix which is chosen in such a way that the second 
order term of the transformed Hamilton function 

h'-^^z) = h'-^^M-^z) (A.ll) 

assumes the particularly simple form 

(z) = Agipi + X:Y(Qi+pi). (A.12) 

k=2 
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Section 2.3 of Ref. [30] provides an explicit procedure for constructing the transfor- 
mation matrix M. 

For the first two steps in the sequence (|A.ip . we actuahy did not give exphcit 
expressions for the generating functions Wi and W2- For conceptual reasons (and to 
justify the notation) we mention that such expressions can be determined (see [30] ) 
but since it is not necessary for the computation we do not discuss these generating 
functions here. The situation is different though for the next steps in (jA.ip which 
rely on the explicit computation of the generating functions Wn with n > 3. In order 
to deal with these higher order transformations we introduce the Poisson bracket of 
two functions A(z) and S(z) which for convenience, we write as 

In this notation the arrows indicate whether the partial differentiation acts to the 
left (on A) or to the right (on B). With the Poisson bracket we associate the adjoint 
operator 

adA-.B^ adAB = {A, B} . (A. 14) 

The transformation (|A.2p then leads to a transformation of the Hamilton function 
jj(n~i) jj(n) -^[ii^ n > 3 which in terms of the adjoint operator ()A.14p reads 

oo ^ 

^^"^ = E^[^dH.j'i?("-^^ (A.15) 

Inserting the power series for H^""^ ' (see Equations (K77\\ and (|X8]) ) in (lAT5]) and 
ordering terms one finds for the s*^ order term of the power series of H^"'^: 

^i"^= E i-A-^wI H^;_lll,y (A.16) 

k=0 

where [-J gives the integer part of a number, i.e., the 'floor'-function. 

Using Eq. (jA.16p one can show that the transformation defined by Eq. (jA.lSp 
satisfies the following important properties for n > 3. Firstly, at step n, n > 3, the 
terms of order less than n in the power series of the Hamilton function are unchanged, 
i.e. 

ijW = , for s<n, (A. 17) 

so that, in particular, = Defining 

P = ad^(2) ={Ff ,.}. (A.18) 

we get for the terms of order n. 



(A.19) 
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This is the so-cahed homological equation which will determine the generating func- 

(n) (n) 

tions Wn for n > 3 from requiring "DHn = 0, or equivalently Hn to be in the 
kernel of the restriction of D to W^. In view of Eq. (jA.lOp this condition yields 

^{n-l) _ g j^g^ _ ^^ 20) 



Section 3.4.1 of Ref. [30] provides the explicit procedure of finding the solution of 
Eq. (|A.20p . In the generic situation where the linear frequencies U2, ■ ■ ■ ,u}d in ()A.12p 
are rationally independent, i.e. m2td2 + . . . + fridiOd = implies m2 = . . . = m^^ = 

(n) 

for all integers m2, . . . , m^, it follows that for odd n, Hn = 0, and for even n, 

i?^") G span {/"i J2"V3"^ ---JT l«l = > (A.21) 

where I = qipi and Jfc = (g^ + p|)/2, with k = 2, . . . ,d. 

Applying the transformation (jA.lSp . with the generating function defined by 
()A.20p . for n = 3, . . . ,N, and truncating the resulting power series at order one 
arrives at the A^**^ order classical normal form (CNF) 

N 

H^cZ(^) = Eo + Y,Hi''Hz). (A.22) 

s=2 

We stress that -ff^NF represents an A^**^ order approximation of the original Hamil- 
tonian H obtained from expressing H in terms of the phase space coordinates 

z(^) = (/>(z) := {cI)Wn ° (PwN-i o • • • o (l>W2 ° 4'Wi){z) (A.23) 

(see (|A.ip ). Truncating the powers series expansion of z*^^) at order A^ we get the 
transformation to the normal form coordinates z = (q, p) = T{q,p) which we intro- 
duced in Eq. ()2.ip . The neighborhood C of validity of the CNF is now defined as a 



neighborhood of the equilibrium zq in which the difference between H^^") and H^p, 
i.e. the remainder term consisting of the nonnormalized tail of terms of order greater 
than A^, and also the difference between the normal form coordinates z and their un- 
truncated version z^^^ can be considered to be sufficiently small for the application 
under consideration (see the comments in Sec. 12. 5p . 

By construction the CNF Hamilton function -ff^NF ^ polynomial of order A^/2 
in the functions / and Jk, see ()2.4p . We thus have 



/ = {/,/?Si} = 0, jk = {Jk,Hg,{} = 0, k = 2,...,d, (A.24) 

i.e. /, and J^, k = 2, . . . ,d, are constants of the motion or integrals for Hamilton's 
equation with Hamiltonian function H^-p. As discussed in Sec. [5] it is the conser- 
vation of these integrals which allows one to construct the phase space structures 
governing reaction dynamics in the neighborhood of the saddle equilibrium point 
Zq. We emphasize that the full algorithm to compute H^^p and the correspond- 
ing coordinate transformation is algebraic in nature, and can be implemented on a 
computer. 
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A. 2 The quantum normal form algorithm 

In the quantum mechanical case we start with a Hamilton operator H which we 
assume to be obtained from the Weyl quantization of a classical Hamilton function 
i7(q, p). Like in the previous subsection, q = {qi,q2, ■ ■ ■ ,Qd) and p = {pi,p2, ■ ■ ■ ,Pd) 
denote the canonical coordinates and momenta, respectively, of a Hamiltonian sys- 
tem with d degrees of freedom. For convenience, we again choose atomic units, 
so that q and p are dimensionless. We denote the corresponding operators by 
q = {qi,q2, ■ ■ ■ ,Qd) and p = {pi,p2, ■ ■ ■ ,Pd)- In the coordinate representation their 
components correspond to multiplication by qj and the differential operators pj = 
—ihesd/dqj. Here ^eff is a dimensionless parameter which corresponds to a scaled, 
effective Planck's constant (see the discussion in Sec.[3|). 

Like in the classical case we will now assume that if(q, p;/ieff) (i-e- the classical 
Hamilton function) has an equilibrium point, zq = (qcPo); of saddle-center-. . .- 
center stability type. Similarly to the classical case the idea of the QNF procedure 
is to approximate the original Hamilton operator if by a simpler Hamilton operator 
where the approximation is valid in a 'neighborhood of the equilibrium zq'. In the 
classical case this simplification is achieved by a symplectic transformation leading 
to a specific choice of the phase space coordinates. In the quantum mechanical case 
the simplification will be achieved by a unitary transformation to a suitable choice 
of a basis for the eigenstates. The basis for computing such unitary transformation 
is the Wigner-Weyl symbol calculus. The crucial point of the symbol calculus is that 
it allows one to identify an operator with its Weyl symbol which is a (in general 
hef[ dependent) phase space functions. This way the QNF theory provides a system- 
atic procedure to obtain a local approximation, Hq^p, of the Hamiltonian if in a 
phase space neighborhood of the equilibrium point zq in order to facilitate further 
computation of various quantities, such as the CRP, of the reaction system under 
consideration, which only depend on properties of H near the equilibrium zq. From 
a computational point of view the symbol calculus is extremely beneficial since the 
resulting phase space functions can be dealt with in a similar way as in the classical 
case described in the previous subsection. 

The Weyl symbol of an operator H is defined as 

ff(q,p;nefT) = y"dx(q-x/2|F|q + x/2)e*^/''-ff. (A.25) 

The map H i— t- ff (q, p; heg) leading to (|A.25p is also called the Wigner map. It is 
the inverse of the transformation which yields a Hamilton operator H from the Weyl 
quantization, Op[-ff], of a phase space function H (the Weyl map) which, using Dirac 
notation, is given by 

H = Op[H] = II -^^ff(q,p;neff) I dx|q-x/2)e-P'^/^^«(q + x/2|. (A.26) 

Accordingly, ff (q, p; hcs) in (|A.25P agrees with the classical Hamilton function ff (q, p) 
in our case. The argument h^s is introduced for convenience since the Weyl symbol of 
the unitarily transformed Hamilton operator we will deal with below will in general 
explicitly depend on hes- 
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Analogously to the classical case we will now construct a sequence of unitary 
transformations which will simply the Hamilton operator order by order in the neigh- 
borhood of the equilibrium point zq. These unitary transformation will be of the form 

Un = e-^'^" , (A.27) 

where Wn is the Weyl quantization Op [Wn] of a symbol Wn which is a homogenous 
polynomial in the space 



W. 



qm 



span {g^i . . . q'^'pf' . ..p^/Ks^ : \a\ + |/3| + 2j = n} . (A.28) 



The main difference between the symbols Wn and the generating functions defining 
the symplectic transformations in the classical case in (|A.2p is that the Wn in (|A.27p 
are polynomials in the phase space coordinates (q, p) and h where for the definition 
of the order, the power of h counts double. 

Using unitary transformations of the form ()A.27p we get the sequence of Hamilton 
operators 

where 

= C>„^("- !)[>*. (A.30) 
Expressing the sequence ()A.29|) in terms of symbols we get a sequence 

H = ii-(o) ^ //(I) ^ i/(2) ^ //(3) ^ . . . ^ /^(^) (A.31) 

which will be the direct analogue of the classical sequence (jA.4p . As we will see 
below, in terms of the symbols the unitary transformations ()A.27p can be viewed as 
h dependent transformations of the phase space coordinates. 

The first two steps of the transformations (|A.3ip are identical to the classical 
case. The first step serves to shift the equilibrium point zq to the origin according to 

(z; h,s) = H^''^ (z + zo; h,s) , (A.32) 

The higher order transformations then work on power series expansions of the symbol 
of the form 

oo 

//(")(z; h,s) = Eo + Y^ Hi^\z; h^s) , (A.33) 

s=2 

where 

An) 



[z,n,s)- 2^ ai\...adW---f3dlf-^' ■■■Id Pi ■■■Pd^oS e W^^. 

\a\ + \l3\+2j=s 

(A.34) 

For n = 1, the coefficients in ()A.34p are again given by the Taylor expansion of H^^'^ 
about the origin 

\ k ,1 — 1 -^l f 
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where we included powers of h only for completeness (i.e. for the type of starting 
classical Hamilton functions in this paper H^^^ does not have terms with nonzero 
powers of h). Like in the classical case the second step of the sequence of transfor- 
mations (|A.3ip consists of a linear transformation of the phase space coordinates by 
a symplectic 2d x 2d matrix M such that the second order term of the symbol 

(z; h,s) = (M-1z; h,s) (A.36) 



takes the form 



k=2 



2 



(A.37) 



The required matrix M is identical to the classical case. 

The reason for the equality of the classical transformations and the quantum 
transformations to second order is the commutativity of the Weyl quantization and 
affine linear symplectic transformations. For the nonlinear transformations corre- 
sponding to steps n > 3 the commutativity ceases to exist and the symbol calculus 
develops its full power. 

In order to simplify the terms of order n > 3 we have to introduce the notion 
of the Moyal bracket. Given two symbols A{z;hes) and B{z;hef[), corresponding to 
operators A and B respectively, the Moyal bracket is defined as 



{A,B} 



M 



-A sin 



^cff 



2 ^^[dq.dp, 



^1 

dpj dqj 



B. 



(A.38) 



The Moyal bracket gives the Weyl symbol of the operator i[A, B]/hcff, where [•,•] 
denotes the commutator. Like in Eq. (jA.lSp the arrows in ()A.38p indicate whether the 
partial differentiation acts to the left (on A) or to the right (on B). Equation ()A.38p 
implies that for h^s — )• 0, 



{A,B}M = {A,B} + Oih,s' 



(A.39) 



i.e. for /igfj = 0, we recover the classical Poisson bracket. Moreover, if at least 
one of the functions A, B is a second order polynomial in the variables (q, p) then 
{A,B}m = {A,B}. Analogously to the adjoint operator associated with the Poisson 
bracket in the classical case (see Eq. (|A.14p ) we now define a Moyal-adjoint operator 
associated with the Moyal bracket: 

MadA : B ^ MadA^ = {A, B}u ■ (A.40) 

The symbol of ^("^ obtained from a unitary transformation of according to 

Eq. ()A.30P is then given by 

= fc! [Madvyj'' /^^""'^ . (A.41) 

k=0 
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Inserting the power series of H^^ and ordering terms one finds for the s^^ order 
terms of the powers series of H^"^^ 

k=0 

where [-J again denotes the floor-function. Equation ()A.42p looks formally the same 
as Eq. ()A.16p in the classical case. The only difference is the occurrence of the 
Moyal adjoint operator instead of the adjoint operator associated with the Poisson 
bracket. This difference however is significant, as the Moyal bracket will in general 
introduce h dependent terms although the original symbol we started with (i.e. the 
classical Hamilton function) had no h dependence. Still, the Moyal bracket has 
similiar properties as the Poisson bracket. Using Eq. ()A.42p we see analogously the 
classical case that at each step n > 3, the terms of order of order less than n remain 
unchanged: 

= , for s<n, (A.43) 

so that we in particular again have -^2"^ ~ ^"i^ ^^^^ classical case. Defining 



P = Mad„(2) ={i?r,-} (A.44) 
we find for the terms of order n 

This is the quantum homological equation which analogously to the classical case 
is used to determine the symbols Wn that lead to the unitary transformations Un 
defined in ()A.27p . Similar to the classical case we determine the Wn by requiring 

(n) (n) 

VHn = 0, or equivalently Hn to be in the kernel of the restriction of V to Wq^^. 
In view of Eq. ()A.45p this condition yields 

Hi^-^) - VWn G Ker V\y,u^^ . (A.46) 

Section 3.4.1 of Ref. ^30j provides the explicit procedure of finding the solution of 
Eq. (|A.46p . Provided the linear frequencies lj2, ■ ■ ■ ,i^d in (|A.37p are rationally inde- 
pendent it again follows that for odd n, = 0, and for even n, 

Hi^^ espan{l"^J^-^J^K..j;^'^Ks^ ■.\a\+j = n/2} , (A.47) 

where I = qipi and Jk = {q\ + P^)/2, with k = 2, . . . , are the analogues of the 
classical integrals. 

Applying the transformation ()A.4ip . with the symbol defined by Eq. (|A.46p . for 
n = 3, . . . , A^, and truncating the resulting power series ()A.33p at the N^^ order one 
arrives at the Weyl symbol Hq^p corresponding to the A^*'^ order quantum normal 



form (QNF) of the Hamiltonian H, 

N 

H^^^^piz; h,s) = Eo + J2 ^s""^ ^eff ) . (A.48) 



s=2 
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We thus see that the computation of the symbol of the QNF operator is very 
similar to the classical case. The major difference is the use of the Moyal bracket in 
the quantum case which is more complicated than the Poisson bracket in the classical 
case. What remains to be done to get the N^^ order QNF operator is to compute 
the Weyl quantization of the symbol Hq^p{z; hcs): 

(A.49) 



^QNF - '-^P 



Trim 
-"QNF 



The Weyl quantization of the classical integrals / and J^, k = 2, . . . ,d, are 

/ = Op[I] = ^{qp + pq), (A.50) 
Jk = Op[Jk] = ^{ql+pl), k = 2,...,d. (A.51) 

Using Eq. ()A.37p and the linearity of the Weyl quantization we get 

d 

Hf^ = \i + Y,LOkJk ■ (A.52) 

k=2 

Since the higher order terms in ()A.48p are polynomials in / and J^, k = 2, . . . ,d (see 
()A.47P ). we need to know how to quantize powers of I and Jk- As shown in [30] this 
can be accomplished using the recurrence relations 

Op [r+i] = /Op [n - 0) ' n^Op [/"-I] (A.53) 



and 



Op [ Jr^] = J,Op [J^] + 0) ' n^Op [rr'] (A.54) 

for k = 2, . . . ,d. Hence, is a polynomial function of the operators / and 

Jfc of the form Eq. (|3.1|) . The coefficients Hn,cx,j are systematically obtained by the 

QNF procedure to compute the symbol H^^p, as desribed above, and the recurrence 
relations ()A.53p and ()A.54p . So like in the case of the classical normal form the 
full procedure to compute is algebraic in nature, and can be implemented 

on a computer. Our software for computing the quantum normal form as well as 
the classical normal form which is recovered for hes = is publicly available at 



http : / /lacms . maths . bris . ac . uk/publications/ software/ index . html, 



In summary the A^**^ order quantum normal form operator is obtained from 

conjugating the original Hamiltonian H by the unitary transformation 

IJ — Q-iWi/hes^-iW2/hcft . . . ^-iWN/ficB 

where we used the fact that the first two steps in the sequence ()A.3ip can also be 
implemented using suitable generators Wi and W2 (see [30] for more details). It 
is an N^'^ order approximation in the sense that the symbol corresponding to the 
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conjugated operator is truncated at order N where the order is defined according 
to Eq. ()A.28p . Using arguments based on the symbol calculus the remainder term 
consisting of the Weyl quantization of the unnormalized tail H^^"^ — Hq^^ is small 
'in a neighborhood of the equilbrium point zq'. The local approximation given by 
the QNF is ideally suited to compute the cumulative reaction probability which is 
directly related to properties of the Hamilton operator in the neighborhood of zq 
because of the simple structure of H^-p. Since is a polynomial function of / 

and Jfc, fc = 2, . . . , (i, we have 

[i,Hi^^^p]=0, [Jk,Hi^Jp] = 0, k = 2,...,d (A.56) 

which is the quantum analog of the invariance of the classical integrals in Eq. ()A.24p . 
The commutativities (|A.56p imply that the eigenstates of Hq^p are product states of 

the one dimensional eigenstates of / and the Jk ■ As we demonstrated in Sections 14.11 
and [5] the simplicity of the spectral properties of I and the upon which the QNF is 
built can be exploited to give an efficient procedure to compute the cumulative reac- 
tion probability and the Gamov-Siegert resonances associated with the equilibrium 

Zq. 
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